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ABSTRACT 


A  triquadratic  isoparametric  solid  element  is  developed  to  study  the  behavior 
of  thick  isotropic  and  laminated  composite  plates.  The  element  is  a  27  noded  La- 
grangian  element  based  on  three  dimensional  elasticity.  Material  characteristics  are 
accounted  by  either  using  laminate  plate  theory  or  three-dimensional  anisotropic 
theory.  Element  matrices  for  nonlinear  stability  analyses  are  derived  based  on  total 
Lagrangian  formulation. 

Results  are  presented  to  compare  with  analytical  solutions  to  validate  the 
elements  behavior.  The  effects  of  various  integration  schemes  on  the  element  per¬ 
formance  are  presented.  Convergence  studies  for  laminated  composites  for  different 
fiber  orientations  are  provided  to  illustrate  applications.  An  analysis  of  thin  plates 
is  carried  out  and  results  for  thick  plates  are  compared  with  available  higher  order 
plate  theories.  One  row  of  elements  in  the  thickness  directions  gives  satisfactory 
results  for  thick  laminates. 
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I.  INTRODUCTION 


A.  OVERVIEW 

The  finite  element  method  provides  a  general  tool  to  solve  problems  of  contin  ja 
such  as  heat  conduction  and  fluid  flow,  but  it  is  most  widely  used  in  structural 
mechanics.  In  structural  mechanics,  the  methodology  is  applicable  for  static  and 
dynamic  response  of  structures  and  in  predicting  the  elastic  stability  limits. 

The  focus  of  the  present  study  is  to  develop  tools  to  analyze  thick  laminatei^ 
composite  plates  and  validate  the  model  by  comparing  with  known  solutions.  More 
specifically,  the  objective  of  the  present  study  is  to  develop  a  finite  element  for  both 
linear  and  nonlinear  analysis  using  three  dimensional  elasticity  relations. 

By  adopting  such  theory  for  thick  plates,  both  isotropic  and  composite,  the 
solutions  account  for  transverse  shear  stresses.  This  approach  eliminates  the  limi¬ 
tations  imposed  by  classical  plate  theory  based  on  Kirchoff-Love  hypothesis  [Batoz, 
1950]  or  higher  order  shear  deformation  theories  [Reddy,  1984,  Lo  et  al.,  1977). 

B.  LITERATURE  REVIEW 

In  this  section,  some  literature  pertaining  to  the  analysis  of  thick  composite 
plates  is  reviewed.  The  finite  element  method  Lm  been  increasingly  used  as  a 
reseeirch  tool,  as  well  as  a  design  analysis  tool,  and  the  methodology  is  rapidly 
evolving  along  with  the  development  of  faster  and  more  efficient  computers.  Basic 
concepts  of  the  theory  of  finite  element  analysis  are  well  documented  [Cook,  et  ru., 
1989].  Yang  (1987)  describes  various  two  dimensional  higher  order  elements  as  well 
as  three  dimensional  solid  elements.  Bathe  (1982)  discusses  the  general  formulation 
of  finite  elements  in  nonlinear  analysis  for  one,  two  and  three  dimensional  elements. 
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based  on  the  total  Lagrangian  formulation  and  the  principle  of  virtual  displacements. 
A  good  source  for  continuum  formulation  may  be  found  in  Malvern  (1969). 

Tsai  and  Pagano  (1968)  establish  a  notation  in  which  composite  lamina  prop¬ 
erties  are  invariant  with  respect  to  lamina  direction.  The  laminate  theory  is  well 
documented  by  Vinson  (1987),  where  the  elasticity  solution  for  “structures  com¬ 
posed  of  composite  materials”  is  given  for  various  cases,  such  as  bending  of  thin 
plates.  Based  on  laminate  theory,  Hoskin,  et  al.  (1986)  outline  procedures  involved 
in  manufacturing  composite  components  and  presents  some  of  its  applications. 

A  higher  order  shear  deformation  theory  of  laminated  composite  plates  was 
developed  by  Lo,  et  al.  (1977).  A  higher  order  nonlinear  theory  of  thick  plates  was 
suggested  by  Reddy  (1984a,  1984b,  and  1985)  and  presented  solutions  (Reddy,  1987) 
and  compared  numerical  results  to  Pagano's  (1969)  elasticity  solution  for  the  case  of 
cylindrical  bending.  Other  elasticity  solutions  are  given  by  Timoshenko  (1951  and 
1959)  and  Eisley  (1989)  who  discusses  the  elasticity  solutions.  The  Heterosis  finite 
element  was  suggested  by  Hughes,  et  al.  (1978)  for  thick  and  thin  plate  bending 
problems.  The  effect  of  reduced  integration  in  isoparametric  finite  elements  was 
presented  by  Zienkiewicz,  et  al.  (1971). 

In  recent  years,  much  work  is  concentrated  on  the  analysis  of  buckling  and  post- 
buckling  response  of  laminated  plates  and  shells  using  nonlinear  analysis.  Ramm 
(1982)  applies  degenerate  finite  elements  to  solve  buckling  of  thin  shells.  Arnold, 
et  al.  (1983)  presents  a  theoretical  analysis  procedure  for  prediction  of  buckling 
and  post-buckling  in  laminated  composite  plates  and  compares  the  results  to  exper¬ 
imental  results.  A  combined  numerical  and  experimental  study  of  the  post-buckling 
behavior  of  composite  panel  is  performed  by  Natsiavas,  et  al.  (1987).  Gujbir  et 
al.  (1989)  use  an  eight  noded  biquadratic  element  to  study  the  effects  of  transverse 
shear  on  the  stability  of  laminated  plates.  Some  solution  algorithms  for  nonlinear 
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analysis  of  structures  by  adapting  modified  Newton- Raphson  and  arc-length  meth¬ 
ods  are  given  by  Kolar  et  al  (1985)  and  Ford  et  al  (1987),  In  the  literature  reviewed, 
there  appears  to  be  no  discussion  on  the  higher-order  solid  element  for  the  analysis 
of  thick  laminated  plates. 

This  research  addresses  the  problem  of  using  a  tri-quadratic  displacement  field 
based  finite  element  based  on  three-dimensional  elasticity  equations.  A  total  La- 
grangian  formulation  is  used  to  derive  relevant  element  nonline2Lr  matrices,  and 
numerical  examples  are  included  to  validate  the  linear  portion  of  the  development. 

Analysis  of  typical  examples  include  slender  bars  under  traction  and  bending 
loads,  thin  and  thick  plates  under  bending  loads  and  effects  of  various  integration 
schemes. 

C.  THESIS  OUTLINE 

This  section  provides  an  overview  of  various  chapters  of  the  thesis.  The  total 
Lagrangian  formulation  for  analyzing  structures  composed  of  three-dimensional  el¬ 
ements  is  presented  in  Chapter  II.  Element  matrices  are  derived  for  both  linear  and 
nonlinear  static  analysis  using  the  incremental  load  method.  The  material  charac¬ 
teristics  account  for  both  linear  isotropic  and  anisotropic  behavior.  Formulas  are 
provided  to  obtain  work-equivalent  loads  for  distributed  body  and  surface  forces. 

Chapter  III  addresses  aspects  of  computational  implementation  of  the  problem 
formulated  in  Chapter  fl.  Test  cases,  example  calculations  and  comparison  with 
classical  solutions  and  other  high  order  theories  are  given  in  Chapter  IV.  Finally, 
Chapter  V  summarizes  the  results  and  reflects  some  suggestions  for  future  work. 
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II.  THEORETICAL  FORMULATION 


A.  INTRODUCTION 

In  this  chapter,  using  the  principle  of  virtual  displacements,  the  stiffness  matrix 
will  be  developed  for  static  equilibrium  of  triquadratic  isoparametric  solid  elements. 
In  the  total  formulation  presented,  both  small  and  large  displacements  are  permis¬ 
sible  for  linear  and  nonlinear  structural  analysis.  For  both  cases,  small  strains  and 
linearly  elastic  material  will  be  assumed. 

The  element  is  developed  for  analysis  of  both  isotropic  amd  composite  struc¬ 
tures. 

B.  GENERAL  DERIVATION  OF  FINITE  ELEMENT  EQUILIBRIUM 
EQUATIONS 

The  principle  of  virtual  work  is  invoked  for  the  general  formulation  of  equilib¬ 
rium  [Bathe,  1982;  Cook,  1989].  The  principle  of  virtual  work  states  that  a  body  is 
in  equilibrium,  if  and  only  if,  the  total  virtual  work  done  by  the  internal  forces  is 
equal  to  the  total  virtual  work  done  by  the  external  forces.  That  is, 

(2.1) 

This  principle  is  equivalent  to  the  minimum  total  potential  energy  principle 
[jllp  =  0],  and  holds  at  any  given  time. 

Consider  a  three-dimensional  body  under  arbitrary  loads  as  shown  in  Figure 
2.1.  Using  a  Cartesian  system,  let  the  loads  be  given  by 

(2-2) 
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{/“} = [ff  /,°  /?r 


(2.3) 


(P)  =  [p.  f;  f;)’'  (2.4) 

where  {/*},  {/®}  and  {F*}  are  surface  tractions,  body  forces,  and  concentrated 
applied  forces  respectively. 

The  displacements  of  a  finite  element  in  the  body  due  to  external  load  is 
denoted  by  {d},  where 

{</}  =  [««  wf  (2.5) 

amd  the  corresponding  strains  are  given  by, 

~  (^*X  Cy,  txy  ]  (^*6) 

for  which  the  corresponding  stresses  are, 

{<r}  =  (or,,  CTyy  <T„  (Ty,  ff,,  (T^]^  (2.7) 

The  total  internal  virtual  work  for  a  finite  element  in  the  body  is  {Se}^ {<7}dv  and 
for  the  whole  body, 

=  J^{S€)^{iT)dv  (2.8) 

where  the  virtual  strains,  {^e},  are 

=  [5c,,  6e„  6t„  6ty,  5e,,  5e,y]^  (2.9) 

The  total  external  virtual  work  is  given  by: 

=  jydY{f^)dv  +  jyd^Y{nds  +  {5«f  }^{F}  (2.10) 
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where  {W*}  denotes  a  surface  displacements  and  represents  point  (nodal) 

displacements  corresponding  to  the  applied  loads,  and  the  virtual  displacements 
{6d}  are 

{6d}^  =  [^u  6ti;]  (2.11) 

On  substituting  equations  2.8  and  2.10  into  2.1,  we  get, 

l^{6t}^{<T}dv  =  J^{6df{f^}dv  +  J^{6d*  finds  +  (2.12) 

It  may  be  noted  that  the  principle  of  complementary  virtual  work  may  have  been 
undertaken,  assuming  small  virtual  stresses  with  true  displacements,  yielding  in  an 
analogous  expression  for  equation  2.12. 

Introducing  the  generalized  Hooke’s  law  for  material  constitutive  relations, 

{<T}  =  [£;]{e}  +  K}  (2.l3) 

where  {<7o}  denotes  the  element  initial  stresses  and  [£]  denotes  the  elasticity  matrix 
of  the  element  material. 

In  general,  the  strain-displacement  relations  are  given  by 

{e}  =  [B]  {d}  (2.14) 

while  the  virtual  strains  are  given  by 

{6€}  =  [B]  {6d)  (2.15) 

Substituting  equations  2.13  and  2.15  into  equation  2.12  and  simplifying,  we  have, 

jm'' {[Bf  IE]IB])  {i}dv  =  l{SdfU‘)dv  +  l{S^f{f}ds 

-  J^{Sdf[B\''{iT.)dv+{Stf{r)  (2.16) 

The  integrations  in  equation  2.16  are  performed  over  the  element  volume  and 
surface,  i.e.,  we  can  evaluate  every  integral  using  the  element  local  coordinates  and 
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assemble  tor  the  global  system  coordinates.  Thus,  we  define  the  global  displacement 
vector  and  the  global  virtual  displacement  vector  as  follows: 


{D}  =  [uiUiWl  U2V2W2  ...UnVnWn]^ 


(2.17) 


(2.18) 


where  n  is  the  total  number  of  nodail  points  in  the  body.  Now  we  define,  for  m 


elements, 


W  =  flliBfmBUv 

11,1  =  JlBf[E\\B]dv 

Jv 


(2.19) 

(2.20) 


where  [A*]  and  [Ar^]  are  the  global  and  local  stiffness  matrices  respectively.  In  addi¬ 


tion,  we  define, 


{«*}  = 

{■■*),  =  liNfU’Uv 

{-•.)>  =  l[Nf{n<^ 

{«()  =  t,f[BfM'B> 

Wh  =  llBf{<r,)dv 


(2.21) 

(2.22) 

(2.23) 

(2.24) 

(2.25) 


(2.26) 


where  {A}  and  {r}  denote  the  global  and  local  load  vectors  and  [A]  and  [A*]  are 
the  displacement  interpolation  (shape  functions)  matrices  for  the  volume  and  surface 
where  traction  is  prescribed.  Using  these  definitions,  we  obtain 

{<£>}'■(*■]  (£>1  =  {«£>)’■({/?.}  + (A) -{«;}  + {r})  (2.27) 


(«)  =  (fla)  +  (fi.) -{«»)  + 


(2.28) 
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By  invoking  the  principle  of  virtual  displacements  and  noting  that  {^£)}  is 
arbitrary,  we  get  the  equilibrium  equations  in  the  following  form: 

[/<]{/)}  =  {/?}  (2.29) 

Equation  2.29  is  the  basic  equation  for  static  equilibrium,  which  also  gives  the 
general  form  for  nonlinear  analysis  with  large  displacements  and  strains. 

C.  INTERPOLATION  SCHEME 

1.  Shape  Functions  (Displacement  Interpolation  Functions) 

In  this  section,  the  interpolation  scheme  for  a  triquadratic  isoparametric 
solid  element  will  be  developed. 

The  one-dimensional  Lagrange  interpolation  function  based  on  parame- 
ters  is  given  by 

P  =  i,  NiPi  =  NiPx  +  iVjPj  +  . . .  +  N,P,  (2.30) 

1X1 

where  Ni,  also  called  the  shape  functions,  are  given  by 

=  (2.31) 

A  triquadratic  solid  element  is  a  three  dimensional  element  in  which  the 
displacements  xi,  v,  and  w  are  interpolated  by  quadratic  langrangian  interpolation 
functions  with  27  nodes.  Figure  2.2  depicts  an  element  in  the  local  non-dimensional 

i 

coordinates  (r,s,t). 

For  an  isoparametric  element,  the  geometry  may  be  interpolated  as, 

2T 

X  =  52 

1-1 

y  =  52 

•xl 
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27  nodes 


(2.32) 


27 

z  =  j]  NiZi 
i=l 

or,  in  a  matrix  form, 


T 

=  (iVKxiyiZi  . .  -  i27y2r^2rl 


where  the  shape  function  matrix  is  given  by 


(2.33) 


'  Ni  0  0  Nj  Q  0  ...  Nrr  0  O' 

[N]=  0  Ni  0  0  N2  0  ...  0  N„  0  (2.34) 

[  0  0  0  0  iVj  . . .  0  0  Nrr  . 

The  shape  functions  and  their  derivatives  in  local  coordinates  are  pre¬ 
sented  in  Appendix  A. 

2.  Jacobian  Transformation  Matrix 

To  obtain  equilibrium  equations  in  the  global  coordinate  system  and  con¬ 
struct  stiffness  matrices,  we  need  the  derivatives  of  the  shape  functions  in  cartesian 
system.  Using  the  chain  rule  for  differentiation,  we  obtain. 


®.r  y,r 

-  -/Vi.,  '  =  -  Ni,y  '  (2.35) 

,  ^i.t  y.»  J  I  ^1.*  . 

where  [J],  the  Jacobian  matrix  is  given  by, 

’  y^r  *.r  ' 

[J]  =  X.,  y.,  z^  (2.36) 

'  *.f  y.i  . 

A  comma  denotes  differentiation,  where  for  example,  iVi.,  =  ^  etc.  Using  the 
shape  function  derivatives,  the  elements  of  the  Jacobian  matrix  may  be  calculated 
and  is  given  in  Appendix  B.  The  global  cartesian  derivatives  may  now  be  obtained 
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where  the  inverse  of  the  Jacobian, 


[r]  =  [J]-‘  (2.38) 

is  given  explicitly  in  Appendix  B. 

D.  STRAIN  DISPLACEMENT  RELATIONS  -  [B] 

1.  Basic  Formulation 

In  this  section,  the  basic  formulation  for  nonlinear  analysis  of  a  general 
solid  body  is  presented  [Refs.  Bathe  (1982),  and  Malvern  (1969)].  First,  some  defi¬ 
nitions  and  notations  will  be  introduced  concerning  the  coordinate  system,  displace¬ 
ment,  stress  and  strain  measures  and  later  on  the  linearized  equilibrium  equation 
will  be  developed  based  on  section  II/B. 

Consider  the  motion  of  a  body,  or  an  element  within,  in  a  fixed  cartesian 
coordinate  system  as  shown  in  Figure  2.3.  We  have  the  body  at  time  0,  t  and  t  +  At 
for  which  the  upper  left  superscript  corresponds.  The  displacements  at  time  t  and 
t  At  are  given  as 

=•  Xi  X,  (2.39) 


Xi  -°x< 

so  that  the  incremental  displacements  are 

=‘+^‘  Ui  -‘«i 

where, 

Ux  =  U  Uj  =  V  U3  =  w 
X^  —  f*  Xj  S  X3  "  t 


(2.40) 


(2.41) 


(2.42) 
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we  use  the  following  notation  for  derivatives  at  time,  say  t  +  At,  with  respect  to 
coordinate  at  time  0  as, 

(2.43) 


= 

0  “«J 


d^Xj 


In  the  present  approach,  we  use  the  total  Lagrangian  formulation,  refer¬ 
encing  all  variables  to  the  undeformed  configuration  at  time  0,  [Ref.  Bathe,  1982; 
Malvern,  1969].  It  is  assumed  that  at  time  0  and  t  the  equilibrium  configuration  is 
known.  Basically,  equation  2.12  needs  to  be  solved  corresponding  to  time  t  +  At. 
Since  we  assume  large  displacements,  and  nonlinear  constitutive  relations  [equa¬ 
tion  2.14],  equation  2.12  may  be  solved  by  incremental  loawl  methods  [Ref.  Ford  k. 
Stiemer,  1987]. 

On  introducing  the  2nd  Piola-Kirchhoff  stress,  it  may  be  shown  that  the 
2nd  Piola-Kirchhoff  stress  tensor  is  energentically  conjugate  to  the  Green-Lagrange 
strain  tensor  [Ref.  Bathe,  1982;  Malvern,  1969]. 

(2.44) 

where  the  2nd  Piola-Kirchhoflf  stress  at  tin^  t  is  defined  as, 

iT 


(2,45) 


such  that  ['x]i  the  deformation-gradient  tensor  is  a  tranformation  operator  from  the 
coordinates  at  time  0  to  time  t.  Note  that  in  the  equation  2.44,  the  right  hand  side 
represents  internal  virtual  work  at  time  At  over  the  volume  at  that  time  while 
the  left  hand  side  has  the  virtual  work  integrated  over  known  configuration  at  the 
reference  volume. 

Assuming  linear  material  behaviour,  we  may  use  the  linear  stress-strain 
relations  (Generalized  Hooke’s  Law)  for  the  2nd  Piolz^Kirchhoff  stress  tensor. 

[■5]  =  (£]  [■<]  (2.46) 
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where  the  Green- Lagrange  strain  tensor  at  time  t  is  defined  as, 

‘Ci;  =  \  (‘“ij  V  ut,. 


(2.47) 


with  i,  j,  A:  =  1,2,  and  3.  On  subs  i,.uting  =‘  u,  -t-  u,,  we  obtain  at  time 


t  +  At, 


4-  2«*.i«fcj 


which  may  be  written  as. 


‘+^‘c  =‘  c  +c 

C|j  —  -r  Ctj 


where,  ‘cjj  is  defined  earlier  and  the  incremental  strain  Cij  is  given  by 


(2.48) 


(2.49) 


c«j  =  -I-  rfij 


(2.50) 


In  matrix  notation. 


U}  =  +  {»?} 


The  linear  incremental  strain  is  identified  as 


Co  =  ^  +‘  +  «*.»  ‘«*j) 


(2.51) 


(2.52) 


in  which  and  *ui,j  are  the  known  displacement  gradients  at  time  t.  The  non¬ 
linear  incremental  strains,  then,  sire  givtn  by 


Vij  —  2“‘'>’ 


(2.53) 


Rewriting  the  equilibrium  equation  as  stated  in  equation  2.12,  using  the  total  La- 


grangian  approach,  we  have, 


/  =‘+^‘  R 


(2.54) 
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where  the  external  virtual  work,  is  assumed  to  be  deformation  independent. 

Using  the  identity  form  equation  2.44,  we  may  write  the  equilibrium  equation  in  the 
undeformed  configuration  as 


/  R  (2.55) 

JOv 

Noting  that  {‘e}  is  displacement  invariant, 

=  {6c}  (2.56) 

The  2nd  Piola-KirchhofF  stress  at  time  t  +  At  may  be  expressed  as 

=  (■>)  +  {^)  (2.57) 

where  {a}  is  the  incremental  2nd  Piola-KirchhofF  stress.  On  substituting  Equations 
2.56,  2.51,  and  2.57  into  2.55  yields, 

+  l{S,,f{‘,)'‘dv  (2.58) 

The  incremental  stresses  are  expressed  using  equations  2.47  and  2.57  as 

{»)  =  [E]  C^'e}  -  (El  {■<)  (2.59) 

which  in  view  of  Equation  2.50  yields 

{,)  =  IE]  {«}  (2.60) 

Referring  to  Equation  2.51  and  neglecting  the  nonlinear  strain  contribution,  we  get 
the  linearized  approximation  as 

{s}  =  [E]{c}  (2.61) 

and 

{Sef  =  [Scf  (2.62) 
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On  substituting  these  into  Equation  2.58  and  rearranging  yields  the  linearized  in¬ 
cremental  equilibrium  equation, 

[  {&)’■  |£|  (e}V«  +  /  [SriVi's^dv  =«‘^'  R-  I  (2.63) 

yOv  JOv  J'^v 

2.  General  Nonlinear  Discretization 

The  general  nonlinear  finite  element  discretization  for  the  27-noded  ele¬ 
ment  is  presented  based  on  the  total  Lagrangian  formulation  discussed  in  the  pre¬ 
vious  section.  Equation  2.52  for  the  linear  incremental  strain  in  cairtesian  form 
yields 

c„  =  +*  iv^ 

®irv  =  “.V  +*  +*  ^.v 

+‘  +‘  U,,  W,  VD^ 

2e„,  =5  -h  w,y  +*  Uy  +  u  „  ‘u.,  +*  y  „  +  v,  ‘v,  +‘  w,y  u>,,  +  w,y  ‘u>^ 

2e,,  =  u.,  +  -h‘  u.,  -f  *u^  -f*  v.,  +  +*  u;^ 

2ety  =  U,y  +  U,y  +  ‘u,,  -f  *  V,y  +  ‘v,,  W,y  -f  * W ,y 

(2.64) 

which  in  matrix  form  is  given  by 

{e}  =  {euo}  -h  {eti}  (2.65) 

The  first  term  on  the  right  hand  side  is  displacement  independent  while  the  second 
term  is  displacement  dependent  with  the  engineering  strains  {c}  represented  by 

{®}  ”  [®**  2er»  2e*|,]  (2.66) 

We  define  the  incremental  displacement  gradient  in  the  global  coordinates  by 

{Ua}  =  [u^  u,y  Vv  w,y  (2.67) 
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Equations  2.64,  2.65,  2.66,  and  2.67  result  in 


{cLo}  =  [Ato]  {C/c} 
{c£,i}  =  [An]  {Ua} 
so  that  [Alo]  and  [i4£,i]  are  given  by. 


(2.68) 

(2.69) 


[Alo]  = 


1  0  0 
0  0  0 
0  0  0 

0  0  0 
0  0  1 
0  1  0 


0  0  0 
0  1  0 
0  0  0 

0  0  1 
0  0  0 
1  0  0 


0  0  0 
0  0  0 
0  0  1 

0  1  0 
1  0  0 
0  0  0 


(2.70) 


[^tl]  = 


‘u.* 

0 

0 

‘v.x 

0 

0 

‘tl,. 

0 

0 

0 

0 

0 

0 

0 

‘W'.v 

0 

0 

0 

0 

0 

‘v.* 

0 

0 

0 

‘u. 

0 

0 

‘u>,^ 

0 

0 

‘t'.x 

0 

‘“.V 

‘ti.. 

0 

‘v.. 

0 

‘lp,J 

0 

(2.71) 


It  may  be  noted  that  the  values  of  are  known  at  the  new  configuration  at  time 
t  +  At.  With  the  displacements  interpolated  by 

«  =  23 

kml 


V  =  ^NkVk 
km\ 

37 

w  -  ^NkWk 

kml 

the  local  displacement  gradients  are  obtained  from 

37 

kml 

37 

kml 

37 

kml 


(2.72) 


(2.73) 
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and  similar  expressions  may  be  attributed  to  v  and  w.  These  gradients  may  be 
represented  as 


r 

0 

0 

N,.r 

0 

0 

^JT.r 

0 

0 

j^i,. 

0 

0 

N7.. 

0 

0 

^27.1 

0 

0 

0 

0 

Nrt 

0 

0 

f^27,t 

0 

0 

'*.1 

0 

Nl.r 

0 

0 

Afj.r 

0 

0 

N27.T 

0 

V,t  >  = 

0 

Ni., 

0 

0 

/Vi.. 

0 

0 

I^27.t 

0 

V.l 

0 

Ny., 

0 

0 

Ni., 

0 

•  •  0 

S27,l 

0 

Nl.r 

0 

0 

N2.r 

0 

0 

N27,t 

Ni., 

0 

0 

N2.. 

0 

0 

1^27., 

Ni., 

0 

0 

Stj  .  . 

0 

0 

N27., 

tlJT 

or  alternatively, 

{UL}  =  [DH]{d} 

where  the  nodal  displacement  vector  is  given  by, 

{d}  =  («i  viWiU^v^  •  •  •  ttfar] 

and  the  local  incremental  displacements  gradient  are  given  by 

{C/l}  =  [Ur  Ut  Vr  V,  W,r 


V27  I 
UfjT  / 

(2.74) 


(2.75) 


(2.76) 


(2.77) 


As  previously  mentioned,  in  isoparametric  finite  elements,  the  same  inter¬ 
polation  functions  are  used  for  approximating  the  geometry  and  the  displacements. 
Using  these  definitions,  we  may  transform  the  displacement  in  global  and  local  co¬ 
ordinates  by  similar  transformations  used  for  the  geometry.  FVirthermore,  we  can 
define  arbitrarily  the  global  and  local  coordinates  to  coincide  at  time  0  configuration, 
thereby,  the  transformation  from  local  coordinates  at  time  0  to  local  coordinates  at 
any  other  time,  say  t  or  f  -H  At  is  identical  to  the  transformation  that  relates  the 
global  coordinates  to  the  local  coordinates  at  any  configuration.  In  other  words,  we 
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can  use  the  same  Jacobian  matrix  defined  previously  for  ail  configurations.  Writing 
the  relation  in  accordance  to  equations  2.37  and  2.38  we  have, 


{Ug}  =  [U  {«U 


such  that, 


[r3l  = 


f(ri  (01  [0)1 
(0)  (n  [0] 
1(01  (0)  [r]J 


and  substituting  equation  2.75  into  2.78  yields. 


{UG}  =  [mDH]{d} 


(2.78) 


(2.79) 


(2.80) 


The  incremental  strains,  then,  may  be  expressed  in  terms  of  nodal  dis¬ 
placements,  and  substituting  equation  2.77  in  2.68  and  2.69 

{<M}  =  (-4K.l|r,llD/fHJ)  (2.81) 

{€(,.}  =  (At, l(r,l(£lffl{i}  (2.82) 

The  strain  displacement  operator  may  be  identified  as 


lBM]  =  (Aullr,)[DAf) 

(2.83) 

lBtil  =  (At,l[r,)l£lffl 

(2.84) 

such  that. 

(Btl  =  (Bwl  +  |Bt.l 

(2.85) 

and, 

II 

(2.86) 

It  should  be  underscored  that  using  only  the  displacement  independent  strains  in 
equations  2.63  and  2.64  results  in  the  linearized  problem  (same  as  linear  small 
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displacement  -  small  strain  formulation)  with 


IBco]  = 


Kz 

0 

0 

A^2.x 

0 

0 

N27^ 

0 

0 

0 

^l.y 

0 

0 

^2.y 

0 

0 

^27.y 

0 

0 

0 

0 

0 

^27.1 

•  •  0 

0 

^27.y 

0 

Nx, 

0 

^2.^ 

N2.y 

0 

^27,t 

^27. V 

A^i,. 

0 

^l,x 

0 

A^2.x  ■  ■ 

^27.1 

0 

^27,z 

^^l.y 

0 

^2.y 

Nl.r 

0 

^27.y 

^27.x 

0 

(2.87) 

The  displacement  dependent  contribution  to  the  strain-displacement  operator 


is  given  by 


IStil  = 


'•«  +'  K, 

/V,^  +'  /V,^ 

ffu  +‘  Wl- 


■w, 

■«,  /V., 

■«,  /V,,  +• /V„ 


'<»,  +' ». 

*•■.  +‘  Nij, 


'“i  +'  “j  ^1* 

'•, /»«.+•■>  Wl- 
■<•,^1,  V.,i¥u. 


^1T, 

'••.t 

A'n- 

■♦''  ■  f  j 

'■i  +' 

Ntit  ♦'  Hvs 


(2.88) 


To  obtain  the  incremental  nonlinear  streun  contribution,  consider  equa¬ 
tion  2.53,  which  has  the  cartesian  components. 


7x* 

= 

+  »"i) 

Vn 

= 

Vmm 

= 

+  v’  +  u;’ ) 

’7f. 

= 

1.  , 
2 

7rx 

= 

1.  , 
2(«.xW.x  + 

’/xy 

= 

1.  . 

-(u^U.,  +  U.xU  y  +  y) 

(2.89) 
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With  the  variations  at  time  t  given  by, 


Sr,„ 

= 

Sv^  ‘u.. 

■  +  Sw^ 

‘u/^ 

Srfyy 

= 

SV^y  ‘Uy 

+  ^ti;.y 

Sviz 

=  ^Ul., 

'^.z  + 

+  Sw^j 

‘^.z 

2Sr)y, 

=  6Uy' 

‘U,y6u,, 

+  ‘ 

*V_y^V,, 

Sw  y  ‘u;,,  -1- 

*Wy6w_t 

2^i7„ 

=  Su^  ' 

+ 

*u^Su^ 

+ 

+  6w^  ‘u;,,  -H 

^w^Swi 

26Tf^  =  6u^  *u^y  +  *u^6u^y  +  6v^  ‘v  ,  +  +  6w^  ‘u;,„  +  *w^6w j^.90) 

(2.91) 

It  is  worth  noting  that  equations  2.90  are  in  exactly  the  same  form  as  the  displace¬ 
ment  dependent  strains  {cn}  given  in  equations  2.64  and  2.65  with  the  incremen|al 
displacements  derivatives  replaced  by  their  variations.  Thus,  we  may  write, 

{Sr,}  =  [Alx]  {SUg}  (2.92) 

iv)  =  [Alt]  {Ug}  (2.93) 

such  that  the  nonlinear  incremental  strain  variation  vector  is  defined  as, 

26f,y,  26Tfg,  26rf^]  (2.94) 

given  in  2.71.  Observing  relation  2.80  for  {Ug},  we  have, 

{SUG}  =  [mDH]{6d}  (2.95) 

and  substituting  for  the  global  variations  into  the  nonlinear  strains  in  equation  2.90, 
we  get 

(«,)’■  =  («)’■  [DHf  [T/  (Auf  (2.96) 
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Using  the  strain- displacement  relations  defined  thus  far,  we  formulate  the 
lineau^ized  incremental  equilibrium  equation,  given  by  2.63  to  take  the  general  form 
as  stated  in  2.29,  to  give, 

-f  [‘*nL])  (2.97) 

where  in  view  of  equation  2.85  is  seen  to  be 

[•k^]=  ljBi.f[E][Bi.fdv  (2.98) 

which  is  the  linear  stiffness  matrix,  and  i  is  the  iteration  number.  This  includes  the 
displacement  dependent  and  independent  contributions. 

When  small  displacements  are  assumed,  the  reduces  to  standard  lin¬ 
ear  stiffness  matrix,  given  by  /o„[5lo]^[^)  [^Loj^dt;.  In  what  follows,  the  derivation 
of  is  described. 

The  2nd-Piola  stress  vector  at  time  t  is  accumulated  such  that, 

{■,}  =  {‘*''■4)  +  |£1  (<)  (2.99) 

Using  equations  2.51,  2.85,  2.86,  and  2.91,  the  incremental  strains  take  the  form 

{€}  =  (2  [/lu]  +  Mwl)  (Fsl  [DH]  {d}  (2.100) 


or  alternatively. 


{£}  =(21Bt,|  +  (B„)){rf) 


(2.101) 


Noting  that  equations  2.99  and  2.100  are  valid  at  any  time  t,  and  by  using  equation 
2.46,  the  2nd  Piola-Kirchhoff  stress  may  be  written  as 


{'»)  =  [£]  (2[Av,]  +  \Alo\)  |r,)113//l  (i)  (2.102) 


and  the  nonlinear  part  of  equation  2.63  becomes,  using  the  relations  2.95  and  2.101, 

({«)’■  [DH]^  irj’'  l/lul’')  1£1 
(2|/l£,,|  +  (/!«,))  ir,]  [DH\  {dfdv  (2.103) 
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we  define  the  nonlinear  strain  contribution,  [Icnl]j 


lkNL]=  [  |Oflr(r,]’'K,|’'i£:i(2l/iL,)  +  (/itol)|r3llfl//|“rft-  (2.104) 

J^v 

It  may  be  recognized  that  [DH]  is  the  strain-displacement  relation  given  by  equation 
2.74  which  corresponds  to  the  linear  contribution  of  the  stiffness  matrix  given  in 
equation  2.97,  and  the  contribution  of  the  non-linear  strains  results  in  the  2nd 
Piola-Kirchhoff  matrix,  [s],  as 

la)  =  iraf  l/ltil’'  (£|  (2  (/It, I  +  |/lu,l)  IPj)  (2. 105) 

[Fa]  is  given  in  equation  2.79.  It  may  be  shown  that  the  matrix  [s)  takes  the  form, 

l^il  [0]  [0] 

[s]=  [Oj  [s,]  (0)  (2.106) 

L  [0]  [0]  [^,] 

The  expression  for  the  second  term  in  the  right  hand  side  of  equation 
2.63  is  evaluated  in  the  same  manner  as  the  linear  and  nonlinear  parts.  On  using 
equations  2.100  and  2.101,  we  obtain, 

/  =  /  {Sdf{2[BufHBi^f)[E]{2[Bu]^Bi^)]{d}°dv  (2.107) 

It  may  be  seen  that  by  defining 

{'"‘‘"f}  =  jf  j2  [Bu]^  +  [Bu]^)  IF)  (2  |Bul  +  |Bi„)|  {<!)"*  (2.108) 

the  expression  {d}^  represents  the  work  done  by  the  external  loads  at  time 

t  -I-  At.  Noting  that 

=  {6e}  (2.109) 

We  approximate  for  the  second  term  in  the  right  hand  side  of  equation 
2.62  such  that 

/  ^  (2.110) 

JOv  -/“v 
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which  represents  the  internal  virtual  work,  so  that  the  right  hand  side  of  the  equi¬ 
librium  equation  2.63  is  the  difference  between  the  external  and  internal  virtual 
work.  On  substitution  of  relations  2.97,  2.102,  2.103,  2.107,  and  2.109  into  2.63 
and  applying  the  principle  of  virtual  displacement  as  shown  previously,  we  arrive 
at  the  incremental  equilibrium  equation  2.96,  which  may  be  solved  by  Newton- type 
methods.  [Ford  and  Stieman,  1987] 


E.  STRESS-STRAIN  RELATIONS 


In  this  section,  the  stress-strain  relations  for  a  composite  material,  to  be  used  in 
the  three  dimensional  analysis  is  developed.  Two  approaches,  one  based  on  classical 
laminate  theory  and  the  other  based  on  anisotropic  material  constitutive  relations, 
are  presented. 

1.  Classical  and  Higher  Order  Laminate  Theories 

Typical  structures  composed  of  composite  materials  are  built  using  sev¬ 
eral  number  of  laminae,  forming  a  laminate.  Each  lamina  consists  of,  typically, 
uniaxial  fibers  embedded  in  a  matrix,  such  as  epoxy-resin,  forming  a  thin  plate. 
Figure  2.4  shows  principal  material  axes,  labelled  1  and  2  in  directions  parallel  to 
and  normal  to  the  fibers,  respectively.  It  may  be  noted  that  in  each  lamina,  there 
exists  a  state  of  plane  stress,  as  shown  in  Figure  2.4. 

Assuming  elastic  orthotropic  material,  (i.e.,  the  lamina  possesses  a  plane 
of  elastic  symmetry  parallel  to  the  x-y  plane),  the  generalized  Hooke’s  law  may  be 
written  as 


Q\i  Qu 

Qi3 

0 

0 

0  ■ 

Cj 

<72 

Q22 

Q23 

0 

0 

0 

C2 

Q33 

0 

0 

0 

4 

£3 

2Q44 

0 

0 

/ 

(4 

symm 

2Qss 

0 

^5 

.  <^6  . 

2C?66  . 

.  ^  . 

(2.111) 
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Figure  2.4:  Lamina  coordinate  system  (2-D) 
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Fig.  2.4;  Lamina  coordinate  system  (2-D) 


where  the  plane-stress  elastic  constants  are  given  by: 


Qn 

Qn 

Q22 

Q44 


Ex 

1  —  Vx2  1^21 

1  -  1/12  1/21 

Q33=T—^ - 

1  —  1/12  i/21 

<^23.  Qss  =  Gi3,  Qee  =  G12 


Qi3  —  Q23  —  O 


(2.112) 


The  subscripts  1,  2,  and  3  correspond  to  normal  stresses  or  normal  strains 
while  4,  5,  and  6  correspond  to  shear  stresses  or  tensorial  shearing  strains  in  yz^  zx, 
and  xy  planes,  respectively. 

The  stresses  in  the  material  coordinate  axes  are  transformed  to  reference 

coordinate  axes  (i,  y,  z)  by  the  following  equation: 

'  0  0  0  — 2mn 

0  0  0  2mn 

0  0  1  0  0  0 

0  0  0  m  n  0 

0  0  0  — n  m  0 

mn  —mn  0  0  0  (m*  —  n^) 

where  the  direction  cosines  m  and  n  are  given  by  m  =  cos  9  and  n  =  sin  0. 

The  strains  may  be  transformed  in  a  similar  manner.  Introducing  the 

strain  transformation,  along  with  equation  2.110  into  equation  2.112  results  in 

QxX  Q\2  0  0  0  2Qi6  Ci 

Qi2  Q22  0  0  0  2Q26  Cy 

_  0  0  Q33  0  0  0  c, 

0  0  0  2Q44  204$  0  €yz 

0  0  0  204$  20$$  0  txz 

^  .  016  026  0  0  0  20$s  .  .  ^xy 

where. 

On  =  +  2  ((5i2  +  2^66)  +  022*^^ 


(2.114) 

k 


Qu  =  (Qii  +  Q22  —  4^66)  +  Qu  (jn*  +  n‘*'j 

Q22  —  Qiin^  +  2  (Qi2  +  2Q6g)  +  (^22^^ 

Q33  =  Q33 

Q\6  =  Qnm^n  -  Q22mn^  -  {Qu  +  ^Qee)  mn 

O26  =  Qnmn^  -  +  (Q12  +  2Qe6)  mn 

^66  =  {<?ii  +  <^22  —  2Qu)  m}n^  +  {m^  — 

Q44  =  Q«m^  +  Qi^n^ 

O45  =  iQss-Q44)mn 

Qss  =  Q44n^ -h  Qssm^  (2.115) 

The  stress  strain  relations  presented  correspond  to  k‘^  lamina.  Now, 
consider  a  laminate  composed  of  N  laminae  for  which,  each  lamina  has  a  different 
orientation  (d),  with  respect  to  the  laminate  x  and  y  axes.  For  linear  eleistic  plates, 
the  function?.!  form  of  the  displacement  may  be  assumed  to  be 

u{x,  y,  z)  =  Uo{x,y)  + zui{x,y) 
u(i,  y,  z)  =  uo(x,y)  +  zvi(x,y) 

w{x,  y)  =  u;o(x,y)  (2.116) 

and,  the  linear  strains  are  given  by, 

=  (2117) 

so  that, 

Cxx  =  l^.x  +  ^«l,x 

tyy  =  UO.V  + 

tzz  =  0 
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(2.118) 


txz  =  2(^1  +  ^0-=f) 

^xy  =  ■*■  ^o,x)  +  ;;^(^i,y  +  vi,i) 

where  uq,  vq  and  wq  are  the  midplane  displacements,  ui  and  vi  are  related  to  the 
rotations  of  the  normals.  It  may  be  noted  that  the  in-plane  strains, 


^zO  no^xi  ^yO  —  ^zyo  —  2^^^*^  no,z) 


(2.119) 


and  the  curvatures  are  given  by 


Kx  —  ni_xi  ^y  —  ^\,yi  ^zy  —  2^^^'^ 


(2.120) 


We  define  stress  resultants  for  plate/shell  type  structures  in  terms,  of 
stresses  and  shears  (see  Figure  2.5)  as  follows  for  the  layer: 

Nx  =  (Txdz 

h 

Qx  =  j\(^yzdz 

Mx  =  /*  (Jxzdz  (2.121) 

•'-7 

Similar  expressions  are  applicable  for  Ny,  Nxy,  Qy,  My  and  Mxy,  where  h  is  the 
lamina  thickness. 

By  summing  all  laminae  over  the  laminate  thickness  in  the  following  man- 


'  f  =  S  {  /  MA  r^Ll  (2122) 

Nxy  J  r-'  [  exyo  U  «xy  J 

which,  in  matrix  form,  may  be  written  as 


(2.123) 
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Figure  2.5:  Stress  resultants  and  couples  in  a  lamina 


Fig.  2.  5;  Stress  resultants  and 

couples  in  a  lamina 


where, 


'4u  =  i;  (4,). 

k=l 

Bi,  =  E  (4,).  (K  -  ftLi)  (2124) 

k=l 

with  i,j  =  1,2,  and  6  and  the  moment  resultants  are  given  by 


{M}  =  [B]{eo}  +  [D]{K}  (2.125) 


where, 

A,  =  J  E  (4i),  (Aj  -  At.)  .2.126) 

with  i,j  =  1,2,  and  6. 

The  displacement  field,  as  stated  in  equation  2.105  is  linear  in  the  thick¬ 
ness  direction,  resulting  in  constant  shear  stresses.  To  get  better  accuracy,  a  higher 
order  displacement  field  may  be  used  [Reddy,  1984]. 

In  order  to  account  for  the  accurate  shear  distribution,  shape  factors  are 
used  in  computing  shear  energies.  These  factors  are  typically  obtained  by  equating 
the  shear  energies.  The  procedure  is  outlined  for  lineair  2ind  cubic  variation  of 
displacement  fields.  The  shear  energy  due  to  transverse  shear  stresses  is  given  by, 

(2.127) 

where  A  is  the  area  bounded  by  the  lamina  surface  dy.  On  using  Hooke’s  Law,  we 
get 

=  j\  (Gxz  +  Gy,  dz  dA  (2.128) 

Equating  this  to  the  linear  displacement  field  and  simplifying, 

U,  =  {^2  Ia  ^  (^1  ^o,»)*]  <^^4^  h  (2.129) 
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Introducing  a  higher  order  displacement  field  yields  a  more  realistic  stress  distribu¬ 
tion,  but  in  doing  so,  a  shape  factor  is  introduced  to  yield  consistent  shear  energy. 

Assuming  a  displacement  field  in  which  the  displacements  are  expanded 
as  cubic  functions  of  the  thickness  coordinate,  while  the  transverse  displacement  is 
assumed  to  be  constant  through  the  thickness,  yields 

W(x,y,i)  =  *^0{r,v)  ^“l(x,v)  "I"  ^  *^2(x,y)  +  ^  *^3(x,y) 

~  ^0(r,y)  ^^l(x,y)  "b  ^  ^2(x,y)  "b  ^  I’Slx.y) 

*"(x.y.x)  =  «>0(r.y)  (2.130) 

With  Uq,  Vo,  and  Wq  being  the  displacements  of  the  midplane,  the  tensorial 
shearing  strains  are  evaluated  as, 

2Cx*  =  [ui(i^y)  +  2u2(x,y)^  "b  3lt3(x,y)2  +  U>o^j 

2Cy,  =  [ui(x.y)  +  2u2(r.y)Z  +  3u3(x.y)2^  +  u?0,y]  (2.131) 

Using  the  condition  that  the  transverse  shear  stresses  vanish  on  the  plate 
top  and  bottom  surfaces,  we  have, 

<^xx  y,  =  <yy*  (x,  y,  =  0  (2.132) 

or, 

Cx*  y,  =  Cyz  y,  =0  (2.133) 

Substituting  relations  2.132  into  2.130,  we  obtain, 

U2  =  V2  =  0 

«3  =  -^(iyo.x  +  «i) 

Vs  =  -3^(«^o,y +  Ui)  (2.134) 
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The  displacement  field  then  becomes 


u 

V 


Uq  +  zui  -  z^—  {Wo^x  +  Ul) 


Vo  +  ZVi  -  Z^—  (lUo.y  +  Vl) 


(2.135) 


W  =  Wq 


(2.136) 


and  the  shearing  strains  are  given  by 

^  (“1  +  ^o,x)  1  -  4 

l-‘'(j)]  (2137) 

Yielding 

Ut  =  {^2  '*■  “’o.r)*  +  Gyj  (vi  +  ti^o.y)^]  j  ^  |l  —  4  ^  j  dz 

(2.138) 

The  shear  energy  ratio  of  the  two  displacement  field  is  found  to  be  so  that  the 
correction  factor  for  constant  shear  stress  using  cubic  displacements  is 

It  is  clear  from  the  discussion  that  the  clctssical  plate  theory  stiffens  the 
plate  by  not  taking  into  account  the  higher  order  terms.  If  we  use  the  higher 
order  theory,  we  need  to  introduce  a  correction  factor  to  the  shearing  strains  of  the 
magnitude  Using  equation  2.113,  the  transverse  shearing  stresses  for  the 
layer  are  given  by 


^XZk  —  ‘^QsSk^XZ  "b  2Q45jCyz 

”  2Q45|,Cr*  +  (2.139) 

and  the  resultants  are  obtained  using  equation  2.122  as 


Qx  —  2  (/Iss^xz  "b  ^45^yz) 
Qy  ~  2(y445Cxi  “b  7444C^j) 


(2.140) 


Note,  that  equations  2.139  and  2.140  are  applicable  for  any  displacement  field. 
Hence,  for  the  higher  order  theory  presented, 

or, 

=  Vt  s  h  -  -  3^1  (*•  -  *‘->)] 

with  i,j  =  4,  5. 

In  the  present  three  dimensional  solid  element,  for  which  only  three  translational 
degrees  of  freedom  per  node  are  defined,  resultants  are  divided  by  the  corresponding 


thicknesses  to  obtain  the  stress-strain  relations  with 


1411  A.12  0  0  0  2i4i6 

1412  i422  0  0  0  2^426 

0  0  0  0  0  0 

0  0  0  2A44  2A,5  0 

0  0  0  2A45  2v455  0 

Ai6  >426  0  0  0  2Ae6 


1^1  =  1 


(2.143) 


For  the  special  case  of  isotropic  material,  the  material  stiffness  matrix  is  given  by 


A-I-2G  A 
A  A  +  2G 


[£]  = 


A  0  0  0 

A  0  0  0 

A-t-2G  0  0  0 

0  GOO 

0  0  G  0 

0  0  0  G 


(2.144) 


where 


(l  +  W(l-2r) 

2.  Three-dimensional  Anisotropic  Theory 


(2.145) 


As  am  alternative  to  using  laminate  theories  to  obtain  [E]  matrix,  we  may 
use  anisotropic  definition  of  the  laminates. 


=  (Oo) 


(2.146) 
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These  relations  eire  approximated  by  obtaining  the  Qij  in  the  laminate  x-y  axes  by 
suitable  transformations  and  transverse  properties  (thickness  direction)  correspond 
to  the  matrix  characteristics. 

F.  CONSISTENT  LOADS 

In  this  section,  we  consider  the  element  noded  loads  vector,  due  to  applied 
loads.  Using  the  virtual  work  principle,  the  distributed  loads,  such  as  surface  loads 
and  body  forces,  are  converted  into  discrete  loads  applied  at  the  element  nodal 
points.  Discretizing  the  distributed  loads  along  these  lines  are  referred  to  as  con¬ 
sistent  or  work-equivalent  loads.  Consider  a  case  where  a  uniform  distributed  load 
acts  on  a  prescribed  face  of  the  element,  as  seen  in  Figure  2.6.  The  consistent  load 
vector  may  be  written  as 

{'■,}= /lAf-rfm  (2.147) 

For  uniform  distributed  surface  loads,  we  have 

{/'}=p(l)  (2.148) 

It  is  worth  noting  that  the  interpolation  functions  on  a  given  surface,  say  t  =  1 
reduces  to  that  of  a  plane  biquadratic  Lagrangian  tscpaxametric  element  and  are 
presented  in  Table  2.1.  Invoking  symmetry,  we  observe  that  the  forces  at  nodes  1, 
3,  5,  and  7  are  equal,  and  similarly,  forces  at  nodes  2,  4,  6,  and  8  are  equal.  On 
using  the  shape  functions  for  the  t  =  1  surface,  we  have, 

5  5’-«  =  5? 

n  =  /j/,  “  ^’’8“  ^’■2  -  (2.149) 
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TABLE  2.1:  SHAPE  FUNCTIONS  FOR  9  NODED  BIQUADRATIC 
ELEMENT 

N,  =  + 

Nz  =  j(l+'-)(l  +  .)-ijV,-iN,- IjV, 
iVi  =  i(l-r)(l+^)-iAr.-iw.-iAr, 

Nz  =  j(l  +  '-)(l-s')-i% 

N,  =  i(i-r^)(i  +  ^)-lw, 

N>  = 

N,  =  5(l-r’)(l-»)-iAf. 

N,  =  (l.r^)(l-s^) 

Note:  Node  numbering  is  referred  to  Figure  2.2  where  t  =  1,  upper  plane. 

Figure  2.6  gives  consistent  element  nodal  loads  for  a  single  element.  As  a  check,  the 
total  pressure  loading  on  the  surface,  2  x  2  x  p  =  4p,  is  seen  to  be  equal  to  the 
sum  of  all  the  discr?^‘/jed  nod2d  point  forces.  The  procedure  may  be  extended  for 
more  than  one  element  by  summing  loads  at  joint  nodes,  as  illustrated  in  Figure  2.7 
for  four  elements. 

G.  INTEGRATION 

In  this  section,  we  summarize  the  Gauss  method  for  numerical  integration, 
including  a  discussion  on  some  aspects  of  integration  schemes. 


36 


\k 

(D 


(D 

H 


Figure  2.7:  Consistent  loads  on  a  mesh 
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Fig. 


1.  Gauss  Quadrature 

The  nature  of  finite  element  matrices  suggests  the  usage  of  numerical 
quadrature.  Gauss’s  integration  scheme  is  the  most  commonly  used  approach  and 
is  adopted  in  the  present  analysis. 

The  method  enables  exact  evaluation  integrals,  consisting  of  polynomi¬ 
als  of  any  order,  by  using  appropriate  order  of  integration.  In  general,  the  Gauss 
quadrature  for  a  function  <^(r,  s,<),  has  the  form 

■f  =  /  /  /  (f>ir,s,t)  dr  ds  dt  f:iY'Y''^WiWjWk4>ir,3,t)  (2.150) 

J-i  J-i  J-i  ijk 

The  integration  limit  reflects  the  limits  of  non-dimensional  ‘master’  isoparametric 
elements,  while  <f>{r,s,t)  represents  the  stiffness  contribution. 

Figure  2.8  demonstrates  the  application  of  the  method  for  a  two  dimen¬ 
sional  biquadratic  element.  Using  the  weighting  factors  as  given  in  Table  2.3,  the 
element  stiffness  matrix  is  evaluated,  for  example,  by  using  a  3rd  order  integration 
scheme  as  follows: 

(A  ]  =  --  {<i>\  -f  <^3  +  </>7  +  <^9)  +  --  {<f>2  +  ^4  +  <^6  +  <t>B)  +  g  g  <^5  (2.151) 

where 

<t>,  =  /i[5(r,s)]^[£;][B(r,s)]  [  J(r,s)  [  (2.152) 

as  is  evaluated  at  Gauss  point  i  as  shown  in  the  Figure. 

2.  Integration  Scheme 

The  term  “full  integration”  refers  to  an  integration  scheme  which  evalu¬ 
ates  the  integral  exaictly  as  shown  in  the  previous  example.  In  the  same  manner,  a 
lower  order  integration  is  referred  to  as  ‘reduced  integration’. 

In  the  present  analysis,  ‘full  integration’  is  used  to  evaluate  the  stiffness 
matrices.  When  a  crude  mesh  is  used,  a  stiffer  structure  is  obtained.  In  general,  there 
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Figure  2,8:  Gauss  points 
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Fig.  2. 8;  Gauss  points 


TABLE  2.2:  SAMPLING  POINTS  AND  WEIGHTS  FOR  GAUSS 
QUADRATURE  OVER  THE  INTERVAL  -1  to  1 


Order  n 


Location  of  Sampling  Point  Weight  Factor  W, 


1 

2 

3 

4 


0. 

±0.57735  02691  89626  =  ±;^ 
±0.77459  66692  41483  =  ±v^ 

0. 

±0.86113  63115  94053  =  ± 

±0.33998  10435  84856  =  ±  [2^]  ^ 
where  r  =  s/T^ 


2. 

1. 

0.55555  55555  55555  =  | 

0.88888  88888  88888  =  | 

0.34785  48451  37454  =  §  -  ^ 

0.65214  51548  62546  =  |  ±  ^ 
and  (2r  —  1)  is  the  polinom  order 


are  two  ways  to  soften  the  structure.  One  way  is  to  refine  the  mesh  and  another  by 
using  ‘reduced  integration’.  Thus,  by  using  a  ‘reduced  integration’  scheme,  a  faster 
convergence  and  more  cost-effective,  accurate  solution  may  be  obtained.  However, 
the  method  suffers  such  drawbacks  as  mesh  instabilities  or  mechanisms,  resulting  in 
a  singular  element  stiffness  matrix. 

H.  BUCKLING  ANALYSIS 
1.  Introduction 

It  is  well  known  that  thin  columns  or  plates  under  axial  compression  tend 
to  buckle.  Elastic  buckling  occurs  when  the  compressive  stress  is  well  below  the 
material  stress  limit.  A  flat  plate  under  axial  compression  shortens  in  the  direction 
of  the  applied  compressive  loads.  This  shortening  results  in  coupling  between  in¬ 
plane  and  out-of-plane  displacements. 
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As  the  applied  compressive  load  increases,  there  is  a  configuration  at 
which  the  plate  offers  no  more  resistance  to  deform,  resulting  in  a  state  of  neutral 
stability.  The  load  corresponding  to  this  configuration  is  referred  to  as  the  buckling 
load  and  constitutes  a  limit  point  on  the  load  response  curve. 

At  this  critical  value,  the  deflection  becomes  very  sensitive  to  any  change 
in  the  configuration.  For  some  structures,  beyond  the  limit  point,  the  load-displace¬ 
ment  path  may  take  any  of  multiple  paths.  The  point  where  the  plate  can  take  any 
of  the  different  paths  is  called  the  Bifurcation  point  and  is  illustrated  in  Figure  2.9. 
In  analyzing  for  nonhnear  response,  the  incremental  load  method  is  adopted,  which 
may  be  summarized  as  follows;  (a)  the  tangent  stiffness  matrix  is  formed,  and  solved 
for  displacements  for  an  incremental  load.  Keeping  the  stiffness  matrix  constant, 
corrections  to  the  incremental  displacements  are  obtained  in  an  iterative  manner 
until  equilibrium  is  achieved,  (b)  total  displacements  for  this  load  are  obtained, 
(c)  a  new  tangent  stiffness  matrix  is  formed  at  this  new  equilibrium  position  and 
steps  (a)  and  (b)  are  repeated.  This  procedure  is  continued  until  the  desired  load 
is  reached  or  the  critical  buckling  load  is  reached. 

2.  Implementation 

In  this  section,  the  Finite-Element  formulation  for  buckling  will  be  pre¬ 
sented  [Bathe,  (1982),  Kolar,  et  al.,  (1985)].  The  problem  of  instability  can  be 
approached  either  by  looking  at  the  equilibrium  of  the  structure  in  the  deflected 
position  and  transforming  all  quantities  to  the  initial  configuration  or  by  solving  the 
system  in  the  current  configuration.  The  former  approach,  described  earher  as  the 
total  Lagrangian  formulation,  is  adopted  here.  By  performing  an  incremental  load 
analysis,  using  the  nonlinear  formulation  described  eairlier,  we  may  write 

([‘A'l]  +  [Knl])  X{P]  -  (2.163) 
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Fig.  2.  9;  Instability  and 

bifurcation  points 


Figure  2.9:  Instability  and  bifurcation  points 


where  {P}  represents  the  total  load  applied  and  is  a  scalar,  referred  to  as 

load  parameter.  The  value  A  scales  the  incremental  load  and  may  be  treated  as  a 
constant  or  variable  during  iterations.  Buckling  load  is  reached  when  displacements 
become  large  with  no  increase  in  the  incremental  load,  i.e.,  the  global  stiffness  of  the 
structure,  as  illustrated  for  a  single  degree  of  freedom  system  in  Figure  2.9,  becomes 
small  and  [K\  tend  to  be  singular.  Thus,  using  Newton- Raphson  and  modified 
Newton- Raphson  methods,  convergence  difficulties  are  encountered  as  buckling  load 
is  approached.  This  is  overcome  by  using  arc  length  methods,  described  in  the  next 
section,  where  the  load  parameter  is  continuously  updated  to  reflect  the  state  of  the 
structure. 

3.  Constant  Arc  Length  Method  [Kolar  and  Kamel  (1985)] 

When  using  the  Newton  type  iteration  schemes,  the  stiffness  matrix  "oe- 
comes  singular  as  limit  points  are  approached.  In  order  to  obtain  post-buckling 
response,  a  method  to  overcome  this  singularity  is  needed.  This  is  accomplished  by 
treating  the  load  parameter  as  a  variable  and  thus  have  an  adaptive  load  incremen¬ 
tation.  This  approach  differs  from  the  conventional  Newton  type  schemes  where  the 
load  level  is  held  constant  for  all  iterations  at  a  given  load  step.  Symbolically,  at 
load  step  m  and  iteration  t,  equation  2.153  can  be  rewritten  <is  follows 

[K]  =  ("*A  +  A(*^  -I-  AA)  {P}  -  (2.154) 

where  [K]  is  the  tangent  stiffness  matrix  at  load  step  m 

A‘+‘  =  A<‘>  -I-  AA  (2.155) 

The  Arc  Length  Method  (ALM)  may  easily  be  visualized  for  a  single  degree  of 
freedom  «is  shown  in  Figure  2.10.  The  displacements  are  updated  as 

="•  {i}  -I-  -I-  Au  (2.156) 
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such  that  corresponds  to  the  displacement  at  the  (i  +  1)‘*  iteration  of  load 

step  m. 

In  the  constant  ALM,  the  radius  of  the  arc  at  each  load  step  is  constant. 
It  is  clear  from  Figure  2.10  that  the  ALM  is  used  in  conjunction  with  the  modified 
Newton- Raphson  method,  and  may  also  be  implemented  with  NR  iteration  schemes. 
The  method  allows  one  to  obtain  postbuckling  response  but  bifurcation  problems 
require  modification  that  will  seek  out  multiple  paths  after  a  limit  point. 

4.  Convergence  Criterion 

For  a  given  load  step,  the  iterations  on  displacements  are  carried  out  until 
a  pre-set  convergence  is  achieved.  There  are  three  convergence  tests  most  commonly 
used,  (a)  Displacement  Convergence,  (b)  Residual  Force  Convergence,  and  (c)  Strain 
Energy  Convergence.  These  criterion  may  be  summarized  «is  follows: 


{Au‘}^{Au'} 

{Au*}^{Au*} 

< 

CtDISP 

iaVU) 

OR.F. 

{AX'f{P}T{P} 

AAi{Aui}{^} 

< 

OtDISPOtR.F. 

(2.157) 


It  may  be  noted  that  {Au‘}  is  the  incremental  displacement  at  iter¬ 
ation,  {y’}  is  the  residual  force  at  iteration,  AA^{P}  is  the  incremental  load  at 
the  first  iteration,  {Au*}  is  the  incremental  displacement  at  the  first  iteration,  and 
q’s  are  the  prescribed  convergence  parameters,  usually  in  the  order  of  10“^  to  10“^. 

It  is  further  noted  that  the  initial  load  {P}  used  to  start  the  analysis  is 
set  arbitrarily  and  only  the  load  parameter  is  modified  automatically  to  go  from 
zero  load  to  the  desired  load  level. 
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III.  PROGRAM  IMPLEMENTATION 


A.  INTRODUCTION 

This  chapter  presents  certain  aspects  of  computer  implementation  of  the  ele¬ 
ments  matrices  developed  in  the  previous  chapter.  As  mentioned  in  Equation  2.25, 
the  general  problem  to  be  solved  is  given  by 


mw  = 

{r} 

(3.1) 

[K]{D]  = 

{R} 

(3.2) 

where  equations  3.1  and  3.2  are  the  static  equilibrium  equations  for  the  element 
and  structural  assemblage  respectively.  If  the  stiffness  matrix,  [K],  is  independent 
of  displacements,  the  analysis  reduces  to  solving  a  set  of  linear  algebraic  equations. 
In  the  case  of  the  stiffness  matrix  being  displacement  independent,  the  structural 
behavior  is  nonlinear,  and  an  incremental  load  analysis  together  with  a  suitable 
iterative  methods  has  to  be  adopted. 

B.  LINEAR  ANALYSIS 

In  the  case  of  linear  analysis,  we  assume  small  displacements  and  small  strains, 
and  the  resulting  force-displacement  relations  are  solved  only  once.  Using  the  dis¬ 
placement  independent  part  in  equation  2.61  to  get  the  strain-displacement  rela¬ 
tions  [Blo],  as  given  by  equation  2.80,  equation  2.20  is  used  to  form  the  element 
stiffness  matrix.  A  series  of  Fortran  subroutines  was  developed  incorporating  the 
element  stiffness  matrix  for  this  element.  The  material  characteristics  may  be  either 
isotropic,  laminate  theory  definitions,  or  anisotropic  description.  The  subroutines 
are  implemented  in  an  existing  computer  program,  FEMCOM,  which  is  capable  of 
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element  assembly  and  subsequently  calculation  of  the  displacement  solution  for  both 
linear  and  incremental  load  methods. 

C.  NONLINEAR  ANALYSIS 

In  this  research  effort,  geometric  nonlinearities,  namely,  large  displacements 
but  small  strains  are  considered.  Consequently,  the  element  consists  of  a  linear 
displacement  independent  stiffness  matrix,  [A’^o],  and  two  other  contributions.  The 
first  contribution  is  due  to  the  linear  displacement  dependent  stiffness  matrix,  [Ali], 
based  on  as  given  in  equation  2.81.  The  other  contribution  comes  from 

nonlinear  stiffness  matrix,  as  given  by  equation  2.95. 

Note  that  the  stress-strain  matrix,  [A],  may  be  used  both  for  isotropic  as  well 
as  composite  materials  using  relations  2.131  and  2.132. 

D.  SOLUTION  PROCEDURE 

1.  Composite  Material 

In  order  to  generalize  the  procedure  of  implementing  the  solid  element 
with  composite  materials,  the  plate  built  of  solid  elements  may  be  stacked  in  all  three 
directions.  Figure  3.1  shows  such  a  stack,  where  rows  of  elements  are  arranged  in  the 
thickness  direction.  For  each  finite  element,  the  stress-strain  matrix  is  computed  in  a 
subroutine  separately,  though  it  would  be  more  efficient  to  compute  it  for  the  whole 
row  of  elements,  taking  into  account  the  appropriate  layers.  In  assigning  a  certain 
number  of  layers  in  each  row,  a  constraint  to  be  noted  is  that  the  total  number  of 
layers  of  all  rows  match  the  number  of  layers  of  the  structure  being  modeled. 

2.  Linear  Case 

As  mentioned  earlier,  the  matrices  corresponding  to  the  linear  displace¬ 
ment  independent  part  was  coded  into  several  subroutines  and  implemented  into  a 
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general  purpose  finite  element  program,  FEMCOM.  The  program  does  automatic 
element  assembly  and  )delds  solutions  to  prescribed  loads.  The  flow  chart  shown  in 
Figure  3.2  shows  various  steps  that  may  be  summarized  ^ls  follows. 

1.  The  material  properties,  model  geometry,  applied  loads,  integration  scheme, 
boundary  conditions  and  solution  parameters  are  input.  The  material  proper¬ 
ties  needed  for  isotropic  material  are  Young’s  modulus  and  Poisson  ratio.  For 
composites,  data  needed  includes  the  number  of  layers,  rows  of  elements,  fiber 
orientations.  Young’s  moduli,  shear  modulus  in  three  directions,  and  Poisson 
ratio. 

2.  Using  the  shape  function  derivatives,  the  coordinates  transformation  relations, 
Jacobian  and  the  strain  displacement  relations  are  established. 

3.  Using  the  specified  Gauss  quadrature,  the  element  stiffness  matrix  is  formed 
in  global  coordinates. 

4.  The  element  global  stiffness  matrix  is  assembled. 

5.  Using  Gauss  elimination  technique,  the  displacement  vector  is  computed. 

6.  Stresses  may  be  computed  using  equations  2.14  and  2.15. 

3.  Nonlinear  Case 

In  order  to  obtain  nonlinear  response,  either  for  studying  the  extension- 
twist-flexure  coupling  or  nonlinear  buckling  and  post-buckling,  the  analysis  pro¬ 
cedure  is  termed  the  incremental  load  method,  and  a  variation  of  Newton-type 
iteration  is  used.  The  element  formulation,  assembly  and  equation  solving  proceed 
as  before,  except  that  additional  element  stiffness  contributions  have  to  be  taken 
into  account.  The  assembly  and  solution  to  get  displacements  needs  to  be  done  as 
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FIG.  3. 2;  Flow  chart- 

linear  analysis 

Figure  3.2:  Flow  chart  —  linear  analysis 
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frequently  as  the  load  steps  increments  and  iterations  continue,  depending  on  the 
solution  strategy  selected.  A  typical  flow  chart  is  given  in  Figure  3.3. 

For  a  given  load  step,  the  incremental  displacements  are  computed  it¬ 
eratively  until  the  convergence  criterion  is  satisfied.  At  that  point,  equilibrium  is 
achieved  and  new  incremental  load  is  applied  and  a  new  tangent  stiffness  matrix  is 
computed.  The  iterations  continue  until  the  new  equilibrium  position  is  obtained. 
In  the  modified  Newton-Raphson  method,  the  tangent  stiffness  matrix  is  kept  con¬ 
stant  for  all  iterations  for  a  given  load  step,  while,  for  the  Newton-Raphson  method, 
the  stiffness  matrix  for  the  whole  structure  is  formed  at  every  iteration.  By  tracing 
the  load-displacement  path,  critical  points,  characterizing  buckling,  and  stable  and 
unstable  regions  of  post-buckling  equilibrium  states  may  be  identified. 
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FIG.  3.3;  Flow  chart- 

Nonlinear  analysis 


Figure  3.3:  Flow  chart  -  Nonlinear  analysis 
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IV.  NUMERICAL  EXAMPLES 


A.  INTRODUCTION  AND  NOTATIONS 

In  this  chapter,  selected  numerical  examples  are  used  to  evaluate  element 
idiosynchrasies  and  demonstrate  its  application  in  solving  critical  structural  com¬ 
ponents  that  use  thick  composites.  Solutions  obtained  here  are  compared  with 
available  elasticity  solutions  or  other  numberical  solutions. 

1.  Material  Properties 

In  all  the  examples  to  be  discussed,  the  material  characteristics  used  are 
cis  follows.  For  isotropic  materials, 

E  =  30x  10®  psi,  1/  =  0.30 

and  for  composite  materials  used  for  laminated  plates,  layer  properties  are  given  by 

El  =  40  X  10®  psi 
E2  =  10®  psi 
Gi2  =  —  0.6  X  10®  psi 

G23  =  0.5  X  10®psi 
V  =  0.25 

An  eight-layered  symmetric  laminate  configuration  using  this  material  is  selected. 
All  the  dimensions  presented  in  this  chapter  are  in  inches.  In  the  discussion  on 
effects  of  numerical  integration  rules,  L  x  M  x  N  notation  refers  to  the  number  of 
integration  points  in  x,  y,  and  z  directions  respectively. 

B.  COLUMNS  AND  BARS 

Two  simple  cases  have  been  selected  as  part  of  the  element  validation  process. 
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1.  Bars 


A  bar  clamped  at  one  end  and  loaded  at  the  other  end  with  uniformly 
distributed  traction  (Figure  4.1a)  was  studied  and  compared  to  the  theory  of  elas¬ 
ticity.  In  the  numerical  solution,  work-equivalent  loads  were  used.  The  dimensions 
of  the  bar  are  10  in.  x  1  in.  x  1  in.,  and  it  is  isotropic.  The  boundary  conditions 
are  given  by 


Figure  4.2  depicts  the  effects  of  reduced  integration  and  mesh  refinement  on  the 
maximum  deflection.  It  is  obvious  that  when  the  mesh  is  refined  in  the  thickness 
direction,  for  instance,  one  element  in  each  of  i  and  y  direction  and  two  in  z  di¬ 
rection  [1x1x2]  mesh,  provides  a  stiffer  solution  than  for  the  [1x1x1]  mesh. 
It  may  be  noted  that  the  full  (F)  and  reduced  {R)  integration  schemes  converge 
to  about  95%  of  the  classical  solution  (See  Appendix  C).  It  may  be  noted  that  the 
classical  elasticity  solution  does  not  account  for  transverse  shear  stresses.  A  reduced 
integration  in  the  axial  direction  (2x3x3)  gives  the  same  results  as  the  full  inte¬ 
gration.  However,  when  reduced  integration  in  the  thickness  directions  is  performed 
(3  x  2  X  2),  the  solution  converges  slowly.  On  using  (2x2x3)  integration  scheme 
in  the  thickness  direction  for  (12  x  1  x  1)  mesh,  spurious  mode  is  observed.  Table 
4.1  summarizes  the  effects  of  various  integration  schemes  and  mesh  sizes. 

2.  Beams 

The  next  example  considered  is  a  clamped,  cantilever  beam  loaded  at  the 
free  end  by  a  shear  load.  Using  the  clamped  boundary  conditions,  dimensions  and 
material  2is  the  previous  example,  solutions  using  full  and  reduced  (/I)  integration 
are  compared  with  the  elasticity  solution  in  Figure  4.3.  The  comparisons  also  include 
the  solution  obtained  using  eight  noded  first  order  solid  element  of  ‘GIFTS’  software. 
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Umax/Uelasticity 


Figure  4.2:  Clamped  Bar  Under  Uniaxial  Load 


TABLE  4.1:  EFFECTS  OF  REDUCED  INTEGRATION  AND  RE¬ 
FINED  MESH,  ON  THE  MAXIMUM  DEFLECTION  OF  CLAMPED 
ISOTROPIC  CANTILEVER  BAR  UNDER  UNIAXIAL  LOAD 


Mesh 

#  d.o.f. 

Integration 

Rule 

Upfiiniiiiimrifl 

1  =  1  X  1  X  1 

5.7 

F  (3  X  3  X  3) 

97.49 

R  (3  X  3  X  3) 

98.17 

RR  (2  X  2  X  3) 

129.56 

111 

F 

99.43 

R 

99.84 

RR 

115.98 

mBwmm 

99 

F 

90.77 

R 

91.20 

RR 

91.83 

3  =  3  X  1  X  1 

165 

F 

100.40 

R 

100.81 

RR 

112.64 

4  =  4  X  1  X  1 

219 

F 

■  iiSB 

R 

■  B 

RR 

■  B 

327 

F 

101.83 

R 

102.80 

RR 

110.74 

12  =  12  X  1  X  1 

F 

■ujm 

■■ 

R 

■■ 

RR 

11,500.00 

U„ 


^ elast 


I  city 


=  Ur, 


AE 

pi 
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The  reduced  integration  shows  a  better  convergence  than  both  the  full  integration 
and  the  first  order  solid  element.  It  may  be  noted  that  the  eight  noded  solid  element 
converges  more  rapidly  than  the  full  integration  scheme  of  the  present  element  up 
to  about  300  degrees  of  freedom  (d.o.f.). 

It  can  be  seen  from  Table  4.2  that  mesh  refinement  in  the  thickness  di¬ 
rection  results  in  reduced  performance  and  one  element  in  the  thickness  direction 
consistently  yields  good  results. 

In  Figure  4.4,  the  effect  of  transverse  shear  deformation  is  studied  for  a 
12  X  1  X  1  mesh  using  reduced  integration  scheme,  and  compared  to  the  theory  of 
elasticity  solution  (See  Appendix  C).  The  results  are  summarized  in  Table  4.3  versus 
the  aspect  ratio. 

It  is  clear  that  for  thin  beams  where  the  elasticity  solution  is  adequtite, 
the  present  element  gives  stiff  solutions,  whereas  for  thick  beams  (^  <  10),  better 
solutions  are  predicted.  The  reason  for  these  effects  may  be  attributed  to  the  trans¬ 
verse  shear  stresses.  In  the  case  of  thin  bars,  or  beauns,  the  element  aspect  ratio 
is  very  large  and  the  parasitic  shear  strains  appear  at  Gauss  points,  resulting  in  a 
phenomena  called  ‘shear  locking’  [Cook,  1989,  and  Hughes,  1978).  When  the  beam 
is  thick  and  the  aspect  ratio  is  of  the  order  of  1/10,  the  transverse  shear  stresses 
start  to  become  significant,  whereas  in  the  elasticity  solution,  they  are  taken  into 
account  only  to  a  limited  degree  together  with  the  restrictions  of  Saint- Venant’s 
principle. 

C.  CLAMPED  PLATES 

An  isotropic  clamped  plate  of  dimensions  20  in.  x  20  in.  x  1  in.  under 
a  central  concentrated  load  is  shown  in  Figure  4.5a.  This  problem  is  studied  for 
mesh  sensitivity  and  the  effects  of  different  integration  rules.  The  present  solution 
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W=Wmax(3EI)/(PL^  )*  100 


Figure  4.3:  End  Loaded  Beam  Bending 
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TABLE  4.2:  EFFECTS  OF  REDUCED  INTEGRATION  AND 
MESH  CONFIGURATION  ON  THE  MAXIMUM  DEFLECTION  OF 
CLAMPED  ISOTROPIC  CANTILEVER  BEAM  LOADED  AT  THE 
END 


Mesh  * 

#  d.o.f. 

tW  =  Wmax  100 

27  solid 

8  solid  t 

(27  solid) 

F 

R 

F 

1  =  1  X  1  X  1 

57 

5.7 

26.7 

9.4 

3  =  3  X  1  X  1 

165 

37.0 

54.0 

47.2 

4  =  4  X  1  X  1 

219 

56.0 

70.3 

60.5 

9  =  9  X  1  X  1 

489 

91.3 

95.4 

85.8 

12  =  12  X  1  X  1 

651 

95.0 

97.3 

89.8 

3  =  1  X  1  X  3 

147 

5.5 

4  =  2x 1 X  1 

135 

16.9 

6  =  2  X  1  X  3 

273 

17.0 

^ela9ticity  “  101-0 

6  =  3  X  1  X  2 

285 

36.0 

9  =  3  X  1  X  3 

399 

39.1 

F:  Full  integration 

•  3  X  3  X  3  for  27  solid 

•  2  X  2  X  2  for  8  solid 

R:  Reduced  integration 

•  2  X  3  X  3  for  27  solid 

t  8  solid  is  generated  in  “GIFTS”. 

*  Mesh  configuration  for  8  solid  is  twice  of  27  solid  in  each  direction,  i.e.,  2x1x3 
for  27  solid  is  4  x  2  x  6  for  8  solid. 
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TABLE  4.3:  CENTER  DEFLECTION  VS.  ASPECT  RATIO  (i)  OF 
AN  ISOTROPIC  CANTILEVER  CLAMPED  BEAM  LOADED  AT  ONE 
END 


— J — 
h 

w 

2 

124.38 

110.07 

5 

103.90 

99.49 

10 

100.98 

97.07 

50 

100.04 

53.43 

100 

100.01 

17.85 

3E/ 

PP 

■"  (I. \h)  Eh 


We 


i  +  2(i  +  .)f 


10’ 


See  Appendix  C. 
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W=Wmax(3EI)/(PL')*  100 


Fig  4.4;  End  Loaded  Beam  Deflection 


Figure  4.4:  End  Loaded  Beam  Deflection 
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is  compared  with  the  elasticity  solution  (See  Appendix  C  for  details  on  elasticity 
solutions). 

Invoking  symmetry  conditions  of  the  problem,  a  quarter  of  the  plate  is  modeled 
with  the  following  boundary  conditions  imposed: 

“  (®’  ^  ("’  ^  (®’  ^  ^ 
u  ^z,  0,  =  u  ^z,  0  ±  =  u;  ^z,  0,  =  0  (4.3) 

and  the  symmetry  conditions. 

“  (^’  (*’  i’  ^  ® 

The  load  was  taken  as  one  quarter  of  the  total  load.  Figure  4.6  shows  the  comparison 
of  a  mesh  composed  of  elements  arranged  in  one  row  of  elements  (N  x  N  x  l)*vs. 
a  mesh  of  the  type  (2  x  2  x  M),  composed  of  M  rows  of  elements  arranged  in  the 
thickness  direction  with  2x2  elements  in  each  row.  Full  integration  is  employed 
in  the  computations.  Table  4.4  summarizes  the  resultant  deflection  and  mesh  sizes. 
It  is  clea*"  from  this  and  the  previous  examples  that  one  row  of  elements  in  the 
thickness  direction  is  2Miequate  to  predict  the  response  of  the  structures.  Figure 
4.7  presents  the  convergence  characteristics  of  three  integration  schemes.  It  may  be 
noted  that  reduced  integration  (3  x  3  x  2)  in  the  thickness  direction  yields  very  close 
results  to  that  of  the  full  integration  scheme.  Reduced  integration  produces  good 
results  by  compensating  for  the  estimation  of  finite  element  approximation.  The 
in-plane  reduced  integration  scheme  (3x2x2)  shows  divergence  in  the  computed 
response.  It  may  be  mentioned  that  using  one  element  to  model  quarter  plate 
resulted  in  much  higher  deflection  than  expected.  This  implies  that  a  one  element 
model  contains  spurious  modes  and  a  one-element  modeling  of  plate/shell  problem 
should  be  avoided.  On  examining  the  convergence  plot,  with  less  than  600  d.o.f.,  the 
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TABLE  4.4:  MESH  COMPARISON  OF  AN  ALL  EDGES  CLAMPED 
RECTANGULAR  ISOTROPIC  PLATE  UNDER  CENTRAL  LOAD  (ip 
=  1000  lb.) 


Mesh 

#  d.o.f. 

w 

1  =  1  X  1  X  1 

37 

13.8529 

4  =  2  X  2  X  1 

145 

4.2230 

9  =  3  X  3  X  1 

325 

5.0690 

16  =  4  X  4  X  1 

577 

5.6592 

8  =  2  X  2  X  2 

275 

4.0133 

12  =  2  X  2  X  3 

405 

3.9955 

16  =  2  X  2  X  4 

535 

3.9597 

18  =  3  X  3  X  2 

591 

5.0647 

a  a 
2  ^  2 

p  a 


h)  Eh^ 

J- -  100 


evaluated  deflection  is  within  90%  of  the  elasticity  solution.  Table  4.5  summarizes 
the  effects  of  various  integration  schemes  and  mesh  sizes. 


D.  SIMPLY  SUPPORTED  PLATES 

Bending  of  a  simply  supported  rectangular  plate  under  uniformly  distributed 
force  is  presented  herein.  (See  Figure  4.5b)  Both  isotropic  and  laminated  plates  are 
investigated  using  one  quarter  of  the  plate,  as  discussed  previously.  The  results  are 
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Fig.  4.6;  Clamped  Plate  Under  Central  Load 


Figure  4.6:  Clamped  Plate  Under  Central  Load 
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Figure  4.7:  Clamped  Plate  Integration  Rules 


TABLE  4.5:  INTEGRATION  RULES  COMPARISON  OF  AN  ALL 
EDGES  CLAMPED  RECTANGULAR  ISOTROPIC  PLATE  UNDER 
CENTRAL  LOAD  (^P  =  1000) 


Mesh 

#  d.o.f. 

w  for  various  integration  rules 

(3x3x3) 

(3x3x2) 

(2x2x3) 

4  =  2  X  2  X  1 

145 

4.2230 

4.2575 

6.9918 

9  =  3  X  3  X  1 

325 

5.0690 

5.0853 

7.6576 

16  =  4  X  4  X  1 

577 

5.6592 

5.6763 

8.1128 

w 


w  (f,  f,  h)  Eh^ 


100 
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compared  with  Classical  Plate  Theory  (CPT)  as  given  in  Appendix  C  and  higher 
order  shear-deformation  (HSDT)  plate  theories  [Reddy,  1985,  Lo  et  al.,  1977]. 

The  following  boundary  conditions  are  imposed, 

w  ^0,  y,  =  w  ^i,  0,  ±1^  =  0  (4.5) 

and  the  symmetry  conditions  are  as  given  in  equation  4.4.  Reduced  integration 
(3  X  3  X  2)  is  adapted  throughout  all  the  computations  presented  in  this  section. 
The  convergence  characteristics  of  the  element  and  comparison  to  CPT  is  depicted  in 
Figures  4.8  and  4.9  for  both  isotropic  and  laminated  plates.  In  the  isotropic  case,  the 
present  element  shows  convergence  within  90%  of  the  elasticity  solution  for  less  than 
200  d.o.f.  In  the  case  of  laminated  plates,  (Figure  4.9),  the  classical  .'>olution  [Vinson, 
1987]  gives  a  more  flexible  solution  than  the  present  element.  Table  4.6  summarizes 
the  deflections  of  the  isotropic  and  laminated  plates  and  mesh  sizes.  It  may  be 
noted  that  the  classical  solution  uses  laminate  theory,  which  neglects  the  transverse 
shear  stresses,  and  hence  the  contribution  of  these  stresses  is  not  taken  into  account. 
This  assumes  more  significance  for  thick  plates  <  10  to  15).  Furthermore,  when 
the  plate  stiffness  in  the  thickness  direction  is  significantly  lower  than  its  stiffness 
in  the  in-plane  direction  and  when  the  shear  modulus  in  the  thickness  direction 
is  significant,  the  classical  laminate  theory  does  not  predict  the  response  of  the 
structure  accurately.  In  the  present  example,  ^  =  20  and  ^  =  40,  G13  «  G12. 
It  may  be  concluded  that  using  laminate  theory  for  bending  of  thick  plates  yields 
a  nonconservative  estimate  of  deflections  and  special  attention  should  be  given  to 
the  stiffness  ratio  amd  shear  modulus  ratio  in  determining  the  response 

of  such  plates.  In  Figures  4.10  and  4.11,  the  maximum  deflection  is  presented  for 
different  aspect  ratios,  of  the  plate.  Both  isotropic  and  laminated  plates  are 
analyzed  and  compared  to  CPT.  In  addition,  the  solution  of  the  laminated  plate  is 
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also  compared  o  Higher  Order  Shear  Deformation  Theory  [Reddy,  1985],  as  shown 
in  Table  4.7. 

As  in  the  beam  bending  r?  oC,  shear  locking  is  observed  for  thin  isotropic  plates. 
For  thick  plates,  (say,  ^  =  4^,  the  computed  deflections  become  significantly  larger 
than  predicted  by  CPT,  as  expected. 

Examining  Figure  4.11,  it  may  be  deduced  that  even  the  HSDT  [Reddy,  1985] 
underpredicts  the  deflections.  For  thin  laminated  plates,  the  shear  locking  effect  is 
not  as  significant  as  observed  for  isotropic  plates.  This  may  be  attributed  to  the 
fact  the  laminated  plate  has  more  flexible  transverse  material  stiffness  in  bending 

isotropic  plate,  so  that  shear  locking  is  expected  to 
plates. 


than  the  coefficients  than  the 
develop  only  for  thin  isotropic 
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TABLE  4.6:  ALL  EDGES  SIMPLY  SUPPORTED  RECTANGULAR 
PLATE,  UNDER  UNIFORMLY  DISTRIBUTED  LOAD 


Mesh 

— 

#  d.o.f. 

Isotropic 

w 

■B9H 

4  =  2  X  2  X  1 

186 

3.2826 

0.3871 

9  =  3  X  3  X  1 

386 

3.6335 

0.3878 

16  =  4  X  4  X  1 

658 

4.0541 

0.4053 

CPT 

4.4335 

*0.3634 

*  Neglects  G13  and  G23.  See  Appendix  C. 

The  uniformly  distributed  load  is  taken  as  the  total  consistent  load  over  the  area  of 
the  quarter  plate. 


F 


The  maximum  deflection  is  taken  at  upper  surface. 

u;(f,  f,  /i)  ^2/1" 
^  - 
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TABLE  4.7:  CENTER  DEFLECTION  VS.  ASPECT  RATIO  OF  SIM¬ 
PLY  SUPPORTED  RECTANGULAR  PLATE  UNDER  UNIFORMLY 
DISTRIBUTED  LOAD 


a 

h 

Isotropic 

w 

Orthotropic 

w 

Reference* 

w 

4 

9.8275 

5.1324 

1.6340 

10 

4.9581 

0.8221 

0.5904 

20 

4.0541 

0.4053 

0.4336 

100 

1.0902 

0.2406 

0.3769 

CPT 

4.4335 

0.3634 

*Reference:  Reddy,  1985. 


7.3 


W=Wniax(E  h^)/(q;H)*  100 


Fii;.  4.H;  .Siniply  Siipiioilcd  Isoimpii:  I’lalc 


Figure  4.8:  Simply  Supported  Isotropic  Plate 


I’ij;.  4.0;  Simply  Supporicd  Laminated  IMatc 


Figure  4.9:  Simply  Supported  Laminated  Plate 


74 


=Wniax(Eh-')/(qa<)'IOO 


Figure  4.10:  Isotropic  Plate  Deflections  vs.  Aspect  Ratio 


Figure  4.11:  Laminated  Plate  Deflection  vs.  Aspect  Ratio 


V.  CONCLUSIONS  AND  SCOPE  FOR 
FUTURE  RESEARCH 


A.  CONCLUSIONS 

This  study  suggests  a  three-dimensional  higher-order  finite  element  to  be  in¬ 
corporated  in  the  analysis  of  thick  plates  composed  of  both  isotropic  and  laminated 
composites.  By  using  a  tri-quadratic  Lagrangian  twenty  seven  noded  solid  element, 
no  assumptions  on  transverse  shear  strains  are  introduced  in  the  formulation.  The 
formulation,  based  on  the  principle  of  virtual  work,  is  presented  for  both  linear 
and  nonlinear  analysis.  The  material  constitutive  relations  for  linear  isotropic  and 
composite  materials  are  presented.  For  composites,  both  laminate  theory  and  three 
dimensional  anisotropic  adaptations  are  described. 

Several  numerical  examples  using  linear  analysis  are  given  for  bars/beams  and 
plates  using  both  isotropic  and  composite  materials.  Three  dimensional  anisotropic 
relations  are  adapted  for  composites.  The  results  show  that  the  present  element  is 
effective  for  analysis  of  thick  beams  and  plates,  but  exhibits  shear  locking  for  thin 
beam  and  plates. 

Spurious  modes  are  revealed  for  single  element  usage  in  plate  modeling,  as  is 
the  case  for  some  other  finite  elements. 

Reduced  Integration  in  the  thickness  direction  for  beams  and  plates  gives  sat¬ 
isfactory  results.  An  interesting  outcome  is  that  one  element  is  sufficient  to  capture 
transverse  deformation  for  thick  laminated  structures  and  mesh  refinement  in  the 
other  two  directions  yields  convergent  solutions. 
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B.  SCOPE  FOR  FUTURE  RESEARCH 

More  numerical  experiments  need  to  be  performed  to  compare  the  present 
solution  to  closed-form  solutions  [Pagano,  1969]  to  evaluate  the  efficacy  of  this  ele¬ 
ment. 

Implementation  of  buckling  analysis  using  the  nonlinear  element  matrices  pre¬ 
sented  herein  is  another  task  that  may  prove  useful  in  predicting  buckling  response 
of  thick  composite  cylinders  subject  to  external  pressure.  By  incorporating  the  cen¬ 
trifugal  force  in  the  external  virtual  work  done  by  body  forces,  this  element  may  be 
used  in  modeling  rotor  blades. 
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APPENDIX  A 

Shape  Functions  and  Derivatives  for  Solid  Element 


Shape  Functions  for  Solid  Element 

Mid-edge  nodes:  Mid-plane  nodes: 


N2=  ir(H-r)(l-s2)t(l-t-l) 

N4=  i(l_r2)«(l-|-s)t(H-t) 

N6  =  -|r(l-r)(l-s2)l(H-l) 

N8  =  -i(l-r2)5(l-s)t(H.O 
Nio  =  -Tr  (1  -h  r)  s  (1  -  s)  (1  -  t^) 
Ni2  =  Ir  (1 -I- r)  s  (1 -fs)  (1 -1^) 
Ni4  =  -Ir  (1  -  r)  s  (!-»-  s)  (1  -  t^) 
Ni6=  |r(l-r)s(l-s)(l-t2) 

N2o  =  -|r(l-hr)(l-s2)<(i_t) 

N22  =  -|(l-r2)s(H-s)t(l-<) 

N24=  ir(l-r)(l-s2)i  (i_t) 

N2s=  lr(l-r^)8(l-s)t{l-t) 


N9=  ^(l-r2)(l-s2)<(l-H) 

N27=-4(l-r2)(l-s2)t(l-0 

Ni,  =  |r(H-r)(l-s2)(l-t2) 
N,5  =  -|r(l-r)(l-s2)(l-t2) 

Ni3=  §(l-r2)8(l-»-s)(l-0 

Ni7  =  -i(l-r2)s{l-s)(l-<2) 


Center  node: 

Ni8  =  (l-r2)(l-s2)(l-t2) 


Corner  nodes: 

Nx  =  l(l  +  r)(l-s)il  +  t)-^-iN2  +  N8  +  Nxo)-\(Nn  +  I^x7  +  N^)-^Nx8 
Nz  =  lil  +  r)(l  +  s)(l  +  t)-^(N2  +  N^  +  Nn)-\(Nn  +  Nx3  +  N2)-lNxs 

^5  =  g(l-r)(l-(-«)(H-t)-i(Ar,-(.Ar6  +  JV,4)-i(A^i3-l-7Vi5-l-7V9)-iNi8 

N7  =  g  (1  —  f)  (1  -  s)  (1  -K)  -  -  {Ne  +  +  A^is)  —  ^  (^15  +  ^17  +  N9)  -  g^is 

Nig  =  g  (1  +  r)  (1  -  s)  (1  -  t)  -  -  (iV20  +  ^26  +  Nio)  —  -  {Nil  +  Ni7  +  N27)  —  g  A^lS 

^21  =  g  (i  +  f)  (1  -h  s)  (1  -  1)  -  -  (iVjo  -1- ^22  +  ^^12)  -  ^  (^11  +  ^13  +  ^^27)  -  -^18 

N2Z  =  g  (1  -  »■)  (1  +  <)  (1  ~  0  ”  2  ^^22  +  ^24  +  A'^m)  -  ^  (A^:3  + -^15  +  A^27)  -  -./Vis 

Nzi  =  g  (1  —  **)  (1  —  *)  (1  —  0  ”  2  ”  4  ~ 
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Shape  Function  Derivatives  for  Solid  Element  •  r  direction 


Mid-edge  nodes:  Mid-plane  nodes: 


N2.r=  il(l-t-2r)(l-S*)(H-<)  Ng,,  =  -rt(l-s2)(l  +  t) 

N<.r  =  -irst(H-s)(l-K)  N27,r=  rt(l-s2)(l-0 

N6,r  =  -}<(l-2r)(l-s2)(n.i)  1  (1 -t- 2r)  (1  -  s2)  (1  _  t2) 

Ns,,  =  lrs<(l-s)(l-K)  Ni5,r  =  -i(l-2r)(l-s2)(l-t2) 

Nio,r  =  -i*  (1  -I-  2r)  (1  -  s)  (1  - 1^)  N,3.r  =  -rs  (1  -I-  s)  (1  -  i^) 

Ni2,r=  is(H-2r)(l-t-s)(l-<2)  N,7.r=  ra(l-s)(l-t2) 

NM,r  =  -|«(l-2r)(H-s)(l-l’) 

Ni6,r=  is(l-2r)(l-s)(l-t2) 

N2o,r  =  —it  (1  -f-2r)  (1  -  (1  —t)  Center  node: 

N22,r  =  -|rst  (1  s)  (1  -  t) 

N24.r=  |t(l-2r)(l-s2)(l-t)  N,8.r  =  -2r(l-s2)(l-l2) 

N26,r  =  -  jrsl  (1  -  s)  (1  -  t) 


Corner  nodes: 


Ni,r  =  g  (1  -  *)  (1  +  <)  -  5  (A^J.r  +  Ns.r  +  A^lO.r)  "  \  {Nn,r  +  A^ir.r  +  A^S.r)  -  g^lS.r 

■Na.r  =  g  (1  +  «)  (H-  <)  ~  2  ^*2.'’)  ~  4  i^n.r  +  ^la.r  +  ^S.r)  -  g^is.r 

Ns,r  —  (1  +  «)  (1  +  <)  -  2  ~  4  ~ 

N7,t  =  -  g  (I  -  s)  (1  + 1)  -  -  (Ne.r  •¥  Ns,,  -t-  Nie.r)  —  ^  (Nu,,  +  Nir,,  +  Ng,,)  - 

^19, r  =  g  (i  “  «)  (1  -  0  ~  2  (^20, r  +  ^26.r  +  ^10, r)  —  (^11, r  +  ^17, r  +  ^27, r)  -  g^l8.>- 

^21, r  =  g  (1  +  «)  (1  ~  0  ~  2  (^20, r  +  ^22.r  +  ^I2,r)  ~  “  (^11, r  +  ^13, r  +  .Ngr.r)  “  g^lSr 

^23, r  =  —  g  (1  +  *)  (1  ~  0  ~  2  (^22, r  +  N^24.r  +  ^M.r)  —  “  (N^13,r  +  Nis,r  +  ^27, r)  “  g'^18.'’ 

iV25,r  =  -g  (1  -  «)  (I  -  0  ~  2  (^24, r  +  Nge.r  +  ^^16,r)  “  “(Nij,,  -t-  Ni7,r  +  N^27,r)  -  g^l«.'‘ 


79 


Shape  Function  Derivatives  for  Solid  Element  •  s  direction 


Mid-edge  nodes; 


Mid-plane  nodes: 


Nj  j  =  -irsl  (1  r)  (1  -f  t) 

N,..=  il(l-r2)(l-h2s)(H-t) 

Ns,*  =  irfi<  (1  -  r)  (1 -1-t) 

Ns,.  =-J  t(l-r2)(l-2s)(l-l-t) 
Nio,.  =  -ir(H-r)(l-2«)(l-t2) 
Ni2,.  =  |r(l-hr)(l-(-2s)(l-<^) 

Ni4,.  =  -}r(l-r)(H.2s)(l-|2) 

Ni6,.  =  5r(l-r)(l-2s)(l-|2) 

N20,.  =  irst  (1 -l-r)  (1 -1) 

N22,.  =  -iMl-'-^)(l  +  2«)(l-0 
N24,.  =  -5rst(l-r)(l-t) 

N26,r=  |t(l-r2)(l-2s)(l-0 


N9..  =  -«<(! -r2)(H-t) 

N27,.  =  st(l-r2)(l-0 

Nu,.  =  -rs  (1  +  r)  (1 
Ni5,.  =  r«  (1  -  r)  (1 -<*) 

Ni3..  =  i(l-r2)(l-h2s)(l-t2) 

Ni7..  =  -^1 -'•*)(! -2«)(l-<") 


Center  node: 

Nis,.  =  -2s(l-r2)(l-l2) 


Corner  nodes; 


^1,.  =  — ^  (1  +  »■)  (1  +  0  ~  2  “*■  ~  4  ^17,.  +  -Ng,.)  — 

^3,.  =  g  (1  +  •■)  (1  +  0  ~  2  ^12,.)  —  ^  (-Nil,.  +  ^13,.  +  ^9,.)  — 

Ny,  =  g  (1  -  »•)  (1  +  0  ~  ^  +  ^8.*  +  ^u.»)  ~  +  ■Nis,.  +  ■Ng,.)  -  g^is,. 

N7,,  =  (1  -  r)  n  -I-  ()  -  i  (iVe,.  +  ^8,.  +  Nie.,)  -  j  (^15,,  +  ^17,.  +  N9,,)  -  ^A^is,. 

A^19,«  =  -g(l  +  '')(l~0“2  (^20,.  +  ^26,.  +  Nl0,$)  —  ^  (^11,.  +  ^17,.  +  ^27,.)  —  g^l8.» 

N^21,«  =  g  (1  +  ’’)  (1  “  0  “  2  (^20,.  +  -^22,.  +  ^12,.)  ~  ^  (^n,.  +  N^13,.  +  ^27,.)  - 

^23,.  =  g  (1  ~  *■)  (1  ~  0  ~  2  (^22,.  +  ^24,.  +  ^14,.)  “  ^  (^13,.  +  ^15,.  +  ^27,.)  —  g^l8,. 

A^25,.  =  “g  (1  ~  *■)  (1  ~  0  -  ^(^24,.  +  ^26,.  +  ^16,.)  “  ^  (^15,.  +  A^17,.  +  A^27,.)  -  g^l8,. 
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Shape  Function  Derivatives  for  Solid  Element  •  t  direction 


Mid-edge  nodes:  Mid-plane  nodes: 


Nj,.  =  lt(l-|-r)(l-s2)(i  +  20  N9.,= 

N4.,  =  5s(l-r2)(l-Hs)(H-2t)  N27.,  =  -i(l -20 

Ns.,  =  -  Jr  (1  -  r)  (1  -  s2)  (1  +  2t)  Nu.«  =  -rt  (1  +  r)  (1  -  s^) 

Ns.t  =-J«(l -r^)  (1-0  (1  +  20  Ni5.,=  rt(l-r)(l-s2) 

Nio,«  =  irst  (1 -t- r)  (1  -  s)  N13,,  = -«<  (1  -  r*)  (1 -|- «) 

Ni2,«  = -|r«t  (1 -I- r)  (!-»-«)  Ni7.«  =  st  (1  -  r^)  (1  -  «) 

Ni4,i  =  |rst  (1  -  r)  (1 -f-s) 

Nj6.«  =  -|rsl  (1  -  r)  (1  -s) 

N20 1  =  —  Jr  (1  -I-  r)  (1  —  s*)  (1  —  20  Center  node: 

N22;«  =  -J«(l-r0  (1  +  0  (1-20 

N2v=  Jr(l-r)(l-s0(l-20  Ni*,.  = -21(1  -  (1  -  *0 

N26.f=  Js(l-r0  (1-0  (1-21) 


Corner  nodes: 


Ni,t  =  g  (1  +  r)  (1  —  s)  —  -  (N2.1  4-  Ns,*  +  Nio,t)  —  j  (^11, «  +  +  ^9,«)  — 

N3.,  =  i  (1  +  r)  (1  -»-  s)  -  i  (^2,.  +  N^,t  +  ^12.1)  -  \  {Nn,t  +  ^13,«  +  N^.,)  - 

=  g  (1  —  r)  (14-  «)  —  -  (A^4,«  +  f^6,t  +  ^M.#)  —  J  (^13.«  +  ^li,t  +  ^9,1)  - 

^7,«  =  g  (1  -  r)  (1  -  s)  -  -  (JVe.i  +  +  Nie.O  —  +  N9,t)  -  g^i8.‘ 

^19.«  =  — ^  (1  +  0  (1  —  0  —  2  (^20,1  +  ^26.«  +  ^10,1)  —  ^  (^11, «  +  ^17, »  +  ^27, t)  — 

^21.«  =  — g  (1  +  r)  (1  4-  «)  —  -  (^20,1  +  ^22.1  +  ^12, «)  —  ^  (^i:,l  +  ^13.1  +  ^27.l)  “ 

^23, «  =  —  g  (1  —  r)  (1  4-  s)  —  -  (A^23,«  +  ^24,t  +  ^  (^13.»  +  ^lS,t  +  ^27, c)  —  g^l8.‘ 

A^25.«  =  (1  -  r)  (1  -  «)  -  -  (^24.«  +  A^26,«  +  ^I6,|)  -  ^  (^15, «  +  ^17,1  +  ^27, l)  -  g^l8.* 
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APPENDIX  B 

Jacobian  Matrix 


Jacobian  matrix  elements: 


27 

Jll  =  X^T  ~  ^  ^  ^i.r 

«=1 
27 

Jl2  =y,T  =  Vi 

i=\ 

27 

Jl3  —  2,r  =  ^i,r  A' 

«=1 
27 

J2I  =X,s  =  21 
«=1 
27 

J22  —  ~ 

i=l 

27 

J23  =  ~  2Z 

»=1 
27 

J31  =  a:,t  = 

«sl 

27 

J32  =  j/,t  =  2Z-^m  y* 

«sl 

27 

J33  =  =  2Z 

«=1 


Elements  of  the  inverse  Jacobian  matrix: 


Fii  =  —  {J27  J33  —  Jxi  Jyi) 

Fij  =  -j  {J\3  Jyz  —  J12  J'Xi) 

Fia  =  -j  {J\2  J23  —  Ji3  J22) 

F21  =  y  ( J23  J31  “  J21  J33) 
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r22  = 

1 

J 

(Jn 

J33  - 

~  Jl3 

J31) 

r23  = 

1 

J 

{J21 

Jl3  ~ 

~  Jn 

J23) 

r3i  = 

1 

J 

{J21  J32  ' 

-  J22 

J31) 

r32  = 

1 

J 

(J12  J31  " 

~  Jn 

J32) 

r33  = 

1 

J 

{J\\ 

J22  ~ 

-  J21 

J\2) 

Jacobian  matrix  determinant: 

J  =  det  [J]  =  Jii  {J22  J33  ~  J23  J32) 

—  J\2  (J21  J33  —  ^23  J31) 

+  «/l3  (*^21  *^32  ~  *^22  *^31) 
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APPENDIX  C 

Theories 


A.  THEORY  OF  ELASTICITY  SOLUTIONS 
1.  Cantilevered  bar  under  traction 


Umax  —  .  _ 


•  P  =  Total  load 

•  L  =  Bar  length 

•  A  =  Cross  section  area 

•  E  =  Young  modulus 

2.  Cantilevered  Beam  under  end  load 


ZEI  21G 

PL^  3 


Reference:  Timoshenko,  1951. 

B.  CLASSICAL  PLATE  THEORY  (CPT) 

1.  All  edges  clamped  rectangular  isotropic  plate  under  central  load 


D 

Eh^ 

12(1-1/2) 

0.00560  for  v  =  0.3 


Reference:  Timoshenko,  1959. 
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2.  All  edges  simply-supported,  rectangular  plate  under  uniformly 
distributed  load 


W„ 


=  a- 


qa 

~D 


Q  =  0  00406  for  v  =  0.3 


3.  Composite 


1  ^  ^  1 
H^mar  =  9^®  II  H 


m=l,3,5...  n=l,3,5... 


r.nD 


D  =  Di\TTi^  -f"  2  (/5,2  4"  {rufi)  4- 1)22^^ 
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TABLE  C-1 

SAMPLE  COMPOSITE  MATERIAL  DATA 


Table  C-1;  Sample  Composite  Material  Data 


INPUT  DATA  ; 


LAMINA; 

THKNES 

THETA  ; 

El 

E2 

;  VI 2 

G12 

8  ; 

0.12500 

0.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

0.60000E4-06 

7  ; 

0.12500 

45.0  ; 

0.40000E+08  : 

O.lOOOOE+07 

;  0.25 

0.60000E-f06 

6  ; 

0.12500 

-45.0  ; 

0.40000E+08  ; 

0.10000E+07 

;  0.25 

0.60000E+06 

5  ; 

0.12500 

90.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

0.60000E+06 

4  ; 

0.12500 

90.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

0.60000E+06 

3  ; 

0.12500 

-45.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

0.60000E+06 

2  ; 

0.12500 

4  5.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

0.60000E+06 

1  ; 

0.12500 

0.0  ; 

0.40000E+08  ; 

O.lOOOOE+07 

;  0.25 

O.eOOOOE+06 

OUTPUT  DATA  ; 


A( i , j )-HATRIX 

0. 12773 E+08  0 . 5594 OE+07-0 . 81226E+06 
0.55940E+07  0 . 17603E+08-0 . 30995E+07 
-0.81226E+06-0.30995E+07  0.59436E+07 

B( i , j ) -MATRIX 

0 .OOOOOE+00-0. 27344E-01  0. OOOOOE+00 
-0.27344E-01-0.16406E+00  O.OOOOOE+00 
O.OOOOOE+00  O.OOOOOE-tOO-0.62500E-01 

D< i , j ) -MATRIX 

0.20743E+07  0.28699E+06  0.72461E+05 
0.28699E+06  0.81544E+06  0.17998E+06 
0.72461E+05  0.17998E+06  0.31612E+06 


Note;  A,B  and  D  matrices  are  evaluated  , neglecting  Transverse  Shear 
contribution  i.e.  G13-G23-0 


Using  Navier's  solution  with  n-m-200,  i.e.  100  terms  for  each 
direction,  as  given  in  the  above,  we  have. 


Hmax  •  0.052328 
3 

Wmax*E  *h 

_ 2 _  *100  -  0.3634 

W  -  J 

q»a 

Where  q»90  ;  a-20 
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